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Abstract 



Systems under external confinement and constraints often show interesting properties. In this 
thesis, we study some systems under external confinement. We begin by finding out the probability 
distribution of end-to-end separation of a semiflexible polymer within the Worm Like Chain model 
using Monte- Carlo simulations (MC). In the constant extension ensemble, where the two ends of 
the polymer are placed in stiff potential traps created by laser tweezers, the probability distribution 
shows triple maxima indicating a non-monotonic force- extension in some intermediate regime of 
polymer stiffness. Whereas this feature is absent in constant force ensemble. Our study on this 
system revealed the ensemble dependence of physical properties for finite sized systems. Then, we 
fix the orientations at the ends of a polymer and find that the orientation, as well as the ensemble, 
control the statistics and mechanical properties of a semi- flexible polymer. We present an exact 
theory to calculate the partition functions in both the Helmholtz and Gibbs ensembles and this 
theory takes care of the orientations at the ends of the polymer. We find multimodality in 
Helmholtz ensemble as a generic signature of semi- flexibility. 

Secondly, we study Laser Induced Freezing (LIF) where a colloidal liquid is constrained by an 
external laser field periodic in one direction. Using a Kosterlitz- Thouless type renormalization 
group calculation and a restricted MC simulation we calculate the phase diagrams for model 
colloids interacting via Hard Disk, Soft Disk and DLVO potentials. The phase diagrams match 
exactly with the corresponding phase diagrams simulated by other groups, thereby proving that 
LIF is indeed a dislocation mediated transition. 

Lastly, we study the phase behaviors and failure mechanism of a two- dimensional solid confined 
within a hard wall channel using MC and Molecular dynamics simulations. This system fails by 
nucleation of smectic phase within the solid. We have shown that thinner strips are stronger! 
The failure is ductile showing reversible plasticity. Density functional arguments can capture 
some of these features qualitatively. We have used a mean field calculation to find out the phase 
diagram of this system in density- channel width plane. We show that fluctuations in quasi one 
dimension lead to very strange behavior, namely a system that looks solid considering its structure 
factor shows vanishingly small shear modulus like a liquid! We study the impact of this reversible 
failure on transport properties. We find that the heat current in response to tensile strain varies 
differently depending on whether the strain is imposed in the confining direction or the other. 
We propose a simple free volume calculation that captures the strain response of heat current, 
exactly, within small strains. 
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1 Introduction 



Every year during the month of March a family of ragged gypsies would set up 
their tents near the village, and with a great uproar of pipes and kettledrums 
they would display new innovations. - G. G. Mdrquez 



Traditionally, the science of thermodynamics and statistical mechanics were concerned with de- 
termining the properties of materials in the so-called thermodynamic limit [1]. In this limit, 
relevant for most materials and experimental situations, the number of particles N ^ oo and 
all statistical mechanical ensembles e.g. Helmholtz (constant number, volume, strain, etc.) and 
Gibbs (constant chemical potential, pressure, stress, etc.) are rigorously equivalent. Also, in 
this limit, for most systems with short ranged interactions and compact shape, the nature of the 
boundaries and boundary fields have no impact on bulk properties. 

Historically, the need for examining systems far removed from the thermodynamic limit arose 
first with the advent of computers and computer simulation techniques in the early 1950s and 
60s. Early computers were not very powerful and the largest system sizes that could be handled 
were small. The need for extrapolating results obtained from small system sizes used in computer 
simulations for predicting phase boundaries, susceptibilities, critical exponents etc. as measured 
in the laboratory spawned the discipline of finite size scaling [2-4]. The early works of Binder 
[5, 6] and Fisher[7-10] are significant in this context. The emerging ideas of the renormalization 
group [11, 12] had a direct impact on this endeavor and aided immensely our understanding of 
the thermodynamics of small systems. 

Simultaneously, technological breakthroughs in semiconductors, magnetic recording devices and 
experiments on surfaces, thin films and adsorbtion pointed out the importance of boundaries 
and surfaces in determining materials properties [13, 14]. Surface phases, wetting, the physics 
of interfaces and surface phase transformations became topics of intense study using computer 
simulations, renormalization group theory and experiments. In the later part of the last century 
rapid advances in primarily two areas of science and technology were responsible for a further 
spurt of activity in trying to understand the role of finite size in determining the properties 
of materials. Firstly, new techniques like the use of lasers in manipulating objects and novel 
microscopic techniques meant that one could measure properties of systems down to the size 
of a large molecule e.g. biologically important polymers like proteins and DNA [15-17]. Rapid 
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advances in biotechnology made it possible for us to ultimately measure forces involved, say, in 
the replication of DNA or in protein synthesis [18, 19]. Nanotechnology on the other hand [20], 
using techniques like laser trapping, atomic force and scanning tunneling microscopy could finally 
make the study of finite systems useful and imperative for its own sake, rather than a precursor 
to taking the thermodynamic limit. 

It is this general context which provides the backdrop of this thesis. In this thesis we have 
focussed on a number of such systems which are either "small" in the sense of being far from the 
thermodynamic limit and/or are acted upon by external fields which produce severe constraints 
that leads to significant change in their behaviour which may be strongly ensemble dependent. 
We study their statistical and mechanical properties, phase behaviors etc. We also study transport 
properties of one such system as it undergoes structural transformations that are controlled by 
external confinement and strains. The structure of this thesis is as follows. 

In the next two chapters we study the properties of a semiflexible polymer. We work within 
a coarse grained model - the worm like chain model - of such a polymer. This model has been 
successful in predicting mechanical properties obtained from single molecule experiments on real 
biological polymers like, DNA, microtubules, actin filaments etc. In chapter-2, we show that in 
the Helmholtz ensemble these polymers can show a non- monotonic force versus extension which 
is an impossibility if the experiment on this finite chain is made in the Gibbs ensemble or in the 
thermodynamic limit. This behavior is obtained in a certain range of stiffness of a semiflexible 
polymer and gives a qualitative signature of semiflexibility vis-a-vis a flexible polymer. In this 
chapter the non- monotonic force- extension curve was obtained for a polymer whose boundaries 
were free to rotate. In chapter-3, we present an exact theory to calculate the properties of a 
semi- flexible polymer for all possible bending rigidities taking into account the orientations of 
ends of a polymer. Using this theory and simulations we establish that imposing constraints 
in boundary orientations vary the statistics and mechanical properties of the polymer. Thus in 
these two chapters we establish that for semi- flexible polymers, both ensembles and boundary 
orientations leave important impact on the system properties. 

In chapter-3, we study the phenomenon of laser induced freezing and reentrant melting. This 
comes about due to constraining a two dimensional system by imposing an external potential that 
is constant in one direction and periodically modulated in the other. We take the simplest case 
where the periodicity of the external potential is commensurate with spacings between the lattice 
planes of the system. Thus the external constraint permeates the whole system. With change in 
the strength of this potential phase transitions occur between a modulated liquid and a locked 
floating solid phase. We use a constrained Monte- Carlo simulation to obtain input parameters 
which are used in a numerical renormalization group scheme to test a recently proposed defect 
mediated melting theory for this system. 

In the following chapter, we study the impact of a different kind of external potential on a two 
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dimensional solid. We confine the solid in a quasi one dimensional hard channel such that, in the 
direction of confinement the solid is only a few atomic layers thick and study its mechanical and 
phase behavior. Study of this system is important due to the recent interest in nano- technology. 
This confinement is a boundary effect and drastically changes the properties of the solid and 
introduces many layering transitions, which are absent in the thermodynamic limit. We also 
study the anomalous and reversible failure of this system under tensile strain. We present mean 
field calculations to substantiate our simulation data. We end this chapter with a discussion on 
enhanced ficutuations in this system due to the reduced dimensionality. 

In chapter-6, we find out the impact of reversible failure on the heat transport properties of 
this system. We use molecular dynamics simulation to calculate the change in heat current as the 
system undergoes tensile strain. The change in heat current in response to strain imposed in the 
confining direction is very different from the case when the strain is imposed in the perpendicular 
direction. We introduce a free volume like theory to calculate the heat current and obtain exact 
match with simulation results up to small strains. 

In summary, in this thesis we have taken up a set of independent problems which are however 
connected by the theme of finite size effects and the effect of boundaries and constraints on the 
determination of overall behavior. Each of the chapters in the thesis are fairly self-contained and 
more or less independent of each other. 



Introduction 



2 Nonmonotonic Force-extension in Semi-flexible Polymer 



"Things have a life of there own", the gypsy proclaimed with a harsh accent. 
'It's simply a matter of waking up their souls". - G. G. Mdrquez 



The simplest model for describing semiflexible polymers without self-avoidance is the so called 
Worm-Like-Chain (WLC) model [21-23]. In this model the polymer is modeled as a continuous 
curve that can be specified by a d— dimensional {d > 1) vector x{s), s being the distance, 
measured along the length of the curve, from one fixed end. The energy of the WLC model is 
just the bending energy due to curvature and is given by 



where u{s) — dx/ds is the tangent vector and satisfies — 1. The parameter k, specifies the 
stiffness of the chain and is related to the persistence length A defined through {u{s).u{s')) — 
e-l^-^'IA. It can be shown that k ^ {d - l)A/2 . 

The thermodynamic properties of such a chain can be obtained from the free energy which 
can be either the Helmholtz (F) free energy or the Gibbs (G) energy. In the former case one 
considers a polymer whose ends are kept at a fixed distance r [one end fixed at the origin and 
the other end at r = (0, ...0, r)] by an average force (/) = dF{r, L)/dr, while in the latter case 
one fixes the force and the average extension is given by (r) = —dG{f,L)/df. It can be shown 
that in the thermodynamic limit L oo the two ensembles are equivalent and related by the 
usual Legendre transform G = F — fr. For a system with finite L/X, the equivalence of the 
two ensembles is not guaranteed, especially when fluctuations become large. We note that real 
polymers come with a wide range of values of the parameter t = L/\ [e.g. A ~ O.l/im for DNA 
while A ~ l/xm for Actin and their lengths can be varied] and fluctuations in r (or /) can be very 
large. Then the choice of the ensemble depends on the experimental conditions. Experiments 
on stretching polymers are usually performed by fixing one end of the polymer and attaching 
the other end to a bead which is then pulled by various means (magnetic, optical, mechanical, 
etc.). In such experiments one can either fix the force on the bead and measure the average 




(2.1) 
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polymer extension, or, one could constrain the bead's position and look at the average force on 
the polymer. In the former case, the Gibbs free energy is relevant while it is the Helmholtz in 
the second case. This point has been carefully analyzed by Kreuzer and Payne in the context of 
atomic force microscope experiments [24]. 

Theoretically, the constant-force ensemble is easier to treat, and infact an exact numerical 
solution has been obtained [25] (though only for t >> 1). In two extreme limits of small force 
and large force WLC model can be solved to obtain force extension relations. In 3D this relation in 
the small force limit is {z)/L — 2Xf/3kBT and in the high force limit is {z)/L — 1 — s/k^TjAXf . 
Data on force-extension experiments on DNA [26] have been explained using this ensemble [25]. 
The case of constant-extension ensemble turns out to be much harder and no exact solution is 
available. The i — > and i — > oo cases correspond to the solvable limits of the hard rod and 
the Gaussian chain. The small and large t cases have been treated analytically by perturbation 
theory about these two limits [27-29]. Numerical simulations for different values of t have been 
reported by Wilhelm and Frey [30], who have also obtained series expansions valid in the small t 
limit. A mean-field treatment has also recently been reported [31]. 

In this chapter we probe the nature of the transition from the Gaussian to the rigid rod with 
change of stiffness as shown by the form of the Helmholtz free energy of the WLC model (or 
equivalently the distribution of end-to-end distance). Extensive simulations are performed in two 
and three dimensions using the equivalence of the WLC model to a random walk with one-step 
memory. We find the surprising result that, over a range of values of t, the free energy has three 
minima. This is verified in a one-dimensional version of the model which is exactly solvable. 

We first note that the WLC model describes a particle in d— dimensions moving with a constant 
speed (set to unity) and with a random acceleration. It is thus described by the propagator 

Z{x, u, L\x', u', 0) 

where in the numerator only paths x{s), satisfying x{0) — x' , x{L) — x, u{Qi) — v! and u{L) — u 
are considered. It can be shown that the corresponding probability distribution W{x, u, L) satisfies 
the following Fokker-Planck equation [27, 28]: 

|^ + 4.V,M/-i-V|M/ = (2.3) 

where V| is the diffusion operator on the surface of the unit sphere in rf— dimensions. The 
discretized version of this model is the freely rotating chain model (FRC) of semiflexible polymers 
[21]. In the FRC one considers a polymer with segments, each of length h = L/N. Successive 
segments are constrained to be at a fixed angle, 9, with each other. The WLC model is obtained, 
in the limit 6',6 — > 0, N ^ oo keeping A = 2b/9'^ and L = Nb finite. In this chapter we report 



j V[x{s)]e-H/''BT 



(2.2) 
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Figure 2.1: Monte-Carlo data for p{v,t) for the 3— dimensional WLC for values of t = 
10(o), 5, 3.33, 2 and 1(V). The inset is a blowup of curves in the transition region {t = 4, 3.85, 3.7) and 
shows the presence of the two maxima. 

the simulation results of this FRC model. In the next chapter we shall see how WLC can also be 
discretized to a Heisenberg spin model with nearest neighbour coupling. 

Here we will consider the situation where the ends are kept at a fixed separation r [with 
a? at the origin and x = f = (0, ...0,r)] but there is no constraint on u and u' and they 
are taken as uniformly distributed. Thus we will be interested in the distribution P{r,L) = 
{S{x — r)) = J duW{f,u,L): this gives the Helmholtz free energy F{r,L) = —Log[P{r, L)]. 
For the spherically symmetric situation we are considering, P{r,L) is simply related to the radial 
probability distribution S{r,L) through S{r,L) = Car'^^^ P{r, L), Cd being a constant equal 
to the area of the rf— dimensional unit sphere. It may be noted that the WLC Hamiltonian is 
equivalent to spin 0{d) models in one dimension in the limit of the exchange constant J — > oo 
(with Jh — K finite) and all results can be translated into spin language. However, for spin 
systems, the present free energy is not very relevant since it corresponds to putting unnatural 
constraints on the magnetization vector. 

2.1 Numerical Simulations 

The simulations were performed by generating random configurations of the FRC model and com- 
puting the distribution of end-to-end distances. To obtain equivalence with the WLC model the 
appropriate limits were taken. We note that because these simulations do not require equilibra- 
tion, they are much faster than simulations on equivalent spin models and give better statistics. 
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Figure 2.2: Monte-Carlo data for p{v,t) for the 2— dimensional WLC for values of t = 

10(o), 5, 3.33, 2 and 1(V). The inset is a blowup of curves in the transition region [t = 4,3.33,2.86) 
and clearly shows the presence of the two maxima. Note that because of ±v symmetry, we have plotted 
data for positive u values only. For the fits at large and small t see text. 



The number of configurations generated was around 10^ for chains of size N — 10^. We verified 
that increasing N did not change the data significantly. As a check on our numerics we evaluated 
(r^) and (r^). Using Eq. (2.3) and following [32] we can compute these (in all dimensions): 



4:kL 8k^(1 



(d-l)L . 

e 2^ 



d-1 



{d 



QAK\d-l) _dL 128K^{d + 5f (^-dl 6Ak^ L{d'^ - 8d + 7) 



2k -\- 



(d-l)L 
2k 



rf3(rf+l)2 (d- 1)2 (d-l)4(rf+l) 

64K\d^ + 23d'^-7d + l) 6AK^L{d^ + 5d'^ -7d+l) 16k^ L\d^ - 3d + 2) 



Infact it is straightforward to compute all even moments, though it becomes increasingly tedious 
to get the higher moments. Our numerics agrees with the exact results to around 0.1% for (r^) 
and 0.5% for {r^) . 

The function P has the scaling form P(r, L) = -^p^r/L, L/X) and we will focus on determining 

the function p{v, t) ^ In Fig. (2.2) and Fig. (2.1), we show the results of our simulations in two and 

three dimensions. At large values of t there is a single maximum at v — r/L — corresponding 

to a Gaussian distribution while at small t, the maximum is close to the fully extended value of 

V — ±1. The transition is first- order-like: as we decrease t, at some critical value, p develops 

^Wilhelm and Prey [30] have looked at the radial distribution S{r, L) which however misses the interesting 
details of the transition. We note that the relevant distribution here is indeed p{v, t) since this gives the 
Helmholtz free energy 
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Figure 2.3: The exact distribution p{v,t) of the 1— dimensional WLC [ Eq. (2.6)] for different 
values of t (10, 5, 3.33, 2, 1). Even for the most stiff chain considered here {t = 1), the distribution has 
a peak at the centre (in addition to the 5— function peaks at ends) though it looks flat. 

two additional maxima at non-zero values of v. Further decreasing t weakens the maximum at 
V — until it finally disappears and there are just two maxima which correspond to the rigid 
chain. 

For the limiting cases of small and large values of t there are analytic results for the distribution 
function and as can be seen in Fig. (2.2,2.1) our data agrees with them. For large t we find that 
Daniels approximation [27], which is a perturbation about the Gaussian, fits the data quite well. 
In the other limit of small t the series solutions provided in [30] fits our data. For intermediate 
values of t neither of the two forms are able to capture, even qualitatively, the features of the 
free energy. Specifically, we note that all the analytic theories(perturbative, series expansions and 
mean-field) predict a second-order-like transition and do not give triple minima of the free energy 
for any parameter value. 

2.2 Exact Calculation in ID 

It is instructive to study the following one-dimensional version of the WLC which shows the same 
qualitative features (the equivalent spin problem is the Ising model). We consider a N step 
random walk, with step-size b which, with probability e, reverses its direction of motion and with 
1 — e, continues to move in the same direction. 

The appropriate scaling limit is: 6 ^ 0, e ^ 0, ^ oo keeping L = Nb, t = L/\ = 2Ne 
finite. Defining Z±(x,L) as the probability of the walker to be at x with either positive or 
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negative velocity, we have the following Fokker-Planck equation: 



T^— T;7t(^+-^-) (2.5) 



dL ' dx ' 2A 
This can be solved for P(x, L) = Z+ Z_ = j^p{x/L, L/X). We get 



e-*/2 

+ ^[<5(^-l) + <^(^ + l)], (2.6) 

where Iq and Ji are modified Bessel functions. In Fig. 2.3 we have plotted p{v,t) for different 
values of stiffness. We find that the probability distribution has three peaks for all values of 
t. Unlike in 2 and 3 dimensions, the 5— function peaks at v — ±1 (which corresponds to 
fully extended chains) persist at all values of stiffness though their weight decays exponentially. 
Similarly the peak at = is always present. 



2.3 Discussion 



The most interesting result of this chapter is the triple minima seen in the Helmholtz free energy 
of the WLC. Physically, this results from the competing effects of entropy, which tries to pull 
in the polymer and the bending energy, which tries to extend it. This form of the free energy 
leads to a highly counterintuitive force- extension curve, very different from what one obtains 
from the constant force ensemble or from approximate theories. In Fig. (2.4) we show the 
force-extension curve for a two dimensional chain with t — 3.33. We see that there are two 
stable positions for which the force is zero. In the constant-force ensemble, it is easy to show 
that d{r)/df — (r^) — (r)^ and so the force-extension is always monotonic. However, in the 
constant-extension ensemble, there is no analogous result {for finite systems), and consequently 
monotonicity is not guaranteed. 

Most of the recent experiments on stretching DNA have t > 100. The distribution is then 
sharply peaked at zero and one expects the equivalence of different ensembles. Experimentally 
the value oH can be tuned by various means, for example, by changing the length of the polymer 
or the temperature. Polymer-stretching experiments can thus be performed for intermediate 
t values. Since we consider the tangent vectors at the polymer-ends to be unconstrained an 
accurate experimental realization of our set-up would be one in which both ends are attached to 
beads [see Fig. (2.5)]. The beads are put in optical traps and so are free to rotate (this setup 
is identical to the one used in refn. [33]). Making the traps stiff corresponds to working in the 
constant-extension ensemble [24] and one can measure the average force. Our predictions can 
then be experimentally verified. 



2.3 Discussion 
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Figure 2.4: The free energy (dotted line) and the corresponding force-extension curve (solid line) 
for a 2— dimensional chain with t = 3.33. 

We make some estimates on the experimental requirements (for a 3-d polymer with stiffness t — 
3.85). Assume that at one end, the origin, the trap is so stiff that the bead can only rotate. We 
make measurements at the other end. The trap-center is placed at fo = (0, 0, zq) and the mean 
bead displacement Az — {{z—zq)) gives the mean force (/) on the polymer. We then consider the 
problem of the polymer in the presence of a trap potential V — [kt{x'^+y'^)+k{z—zoy]/2. Assume 
kt » k so we can neglect fluctuations in the transverse directions. The distribution of the bead's 
position in the presence of the potential is given by Q{r) = e'l^^^^^+^^^y J dVe-'^l^^^^+^^^^l . 
For a stiff trap, we can expand F about f—fo and find the average displacement of the bead. 

dF 1 d'^F 

F{r) = F(fo) + g^lroixi - Xq,) + 2d^^.\^o'<''' ~ ''o^^^''^' ~ ""o^)" 

Due to spherical symmetry F{r) — F{r). The operator identity ^ — ■§^-§^ — leads to 

axi r or 



d^F 
dxidxj 



IdF ^ 
r or 



2C 'iCC j 


fd'^F 






dF 


— dF 


' dx 


dy 



r dr ) 



d'^F _ a'^F 
dx'^ dy'^ 



Fq/zq = Go and 



d^F 
dx'^ 



F" Therefore, 



d^F 



= for i 7^ j, 



F{r) + V{r) = Fo + F^^z - zo) + \{Go + h){x^ + y^) + \{F^ + k){z - zof 
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« Zo * 

Figure 2.5: A schematic of the experimental set-up required to realize the constant-extension 
ensemble discussed in this chapter. For a stiff trap the average displacement of the bead (Az) from the 
trap center is small and the average force on the polymer is: (/) = —k{Az). 

writing F(ro) = Fq. Then the average displacement can be written as, 

where k' — k + F"{zq). In short, we have shown that = J d^r{z — zo)Q{r) — —{f)/k', 
where k' — k + F"{zq) ^ k (valid except when zq L) and (/) — Fq — F'{zo). The 
rms fluctuation of the bead about the trap center is given by z^^g — ksT/k. Hence we get 
Az = -{f)z',^J{kBT) = -{{f)L/kBT){zfUL). The scaled force {f)L/{kBT) is of order 
0.1. The different minima are separated by distances ^ 0.2L, hence to see the effect we need to 
have Zrms/L < 0.1. Thus finally we find that the typical displacement of the bead Az is about 
O.OlZrms- This is quite small and means that it is necessary to collect data on the bead position 
over long periods of time. 

As suggested in [30], a more direct way of measuring the Helmholtz free energy would be to 
attach marker molecules at the ends of the polymer and determine the distribution of end-to-end 
distances. Fluorescence microscopy as in [15] could be another possible method. It is to be 
remembered of course that real polymers are well-modeled by the WLC model provided we can 
neglect monomer-monomer interactions (steric, electrolytic etc.). Thus the experiments would 
really test the relevance of the WLC model in describing real semiflexible polymers in different 
stiffness regimes. 



2.4 Conclusion 
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2.4 Conclusion 

In conclusion we have presented some new and interesting properties of the WLC model and have 
pointed out that polymer properties are ensemble-dependent. This is a finite size effect. In this 
chapter we have given one example of qualitative differences in force-extension measurements in 
different ensembles. Other quantitative differences will occur even in more flexible chains and 
should be easier to observe experimentally. If the ends of the polymer are free to rotate but 
the end to end separation is constrained to be at a constant distance, such that the system is in 
constant extension ensemble, then the free energy of the polymer is shown to have a triple minima 
structure for some intermediate values of stiffness. In the next chapter we shall investigate the 
robustness of this triple minima structure as some fixed orientations. 
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3 Semiflexibile Polymer: Ensemble and End- Orientation 



The children would remember for the rest of their lives the august solemnity 
with which their father, devastated by his prolonged vigil and by the wrath of 

his imagination, revealed his discovery to them: "The earth is round like an 
orange". - G. G. Mdrquez 



In the last chapter we have studied the statistical and mechanical properties of a semiflexible 
finite polymer with its ends left free to rotate. The bending rigidity coupled with the finite size 
of the polymer gave rise to inequivalence of ensemble and a very interesting triple minima in 
free energy in the Helmholtz ensemble and consequently non- monotonic force extension which is 
absent in Gibbs ensemble. In this chapter, we shall discuss an exact theory that gives results which 
match exactly with the simulations. Further, we shall go over to different boundary conditions by 
fixing the orientations at the ends of a polymer and study its impact, which is very non- trivial, 
on the end to end probability distributions and mechanical properties. 

Microtubules and actin polymers constitute the structure of cytoskeleton that gives shape, 
strength and motility to most of the living cells. They are semiflexible polymers in the sense that 
their persistence lengths A are of the order of their chain lengths L such that t — L/X\s small and 
finite. For example, Actin has A = 16.7 iim, L ~ 30 //m[15, 16], Microtubule has A = 5.2 mm 
and statistical contour length L ~ 10 /im[16], double stranded DNA has A = 50 nm and contour 
length L ~ 300 nm[18]. While it is obvious that in the thermodynamic limit of t — > oo, Gibbs 
(constant force) and Helmholtz (constant extension) ensemble predict identical properties, the 
same is not true for real semi- flexible polymers which are far away from this limit. In biological 
cells actin filaments remain dispersed throughout the cytoplasm with higher concentration in the 
cortex region, just beneath the plasma membrane. Microtubules, on the other hand, have one 
end attached to a microtubule- organizing centre, centrosome in animal cells. Another polymer, 
microtubule- associated proteins (MAP) attach one or both their ends to microtubules to arrange 
them in microtubule bundles [17]. Thus, end point orientation of polymers play a crucial role 
in many important phenomena. For instance, gene- regulation in the cell is controlled by DNA- 
binding proteins, many of which loop DNA with fixed end orientations [34-36]. Thus it becomes 
important to understand the statistics and the mechanical properties of semi- flexible polymers 
with different possibilities of end orientations and ensembles. 
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Figure 3.1: Three possible experimental set- up for force- extension measurements, (a) left and 
right ends are held by optical traps, (b) left end is anchored to a surface and the right end is held by 
optical trap, (c) left end is anchored to a surface and the right end is held by the functionalized tip 
of an AFM cantilever, zq denote the position of the free dielectric bead (for optical trap) or the free 
cantilever tip in absence of the polymer, Az denote the displacement associated with the force exerted 
by the polymer. 

During the last decade many single molecule experiments have been performed on semi- flexible 
polymers[18, 19, 26, 37]. This has been done by using optical tweezers[37] , magnetic tweezers[38], 
AFMs[39] etc. In optical tweezer experiments one end of a polymer is attached to a dielectric 
bead which is, in turn, trapped by the light intensity profile of a laser tweezer. In this case the 
dielectric bead is free to rotate within the optical trap. On the other hand, attaching an end of a 
polymer to a super-paramagnetic bead, one can use magnetic field gradients to trap the polymer 
using a magnetic tweezer setup. In this case one can rotate the bead while holding it fixed in 
position by changing the direction of the external magnetic field. In AFM experiments on the 
other hand one end of polymer is trapped by a functionalized tip of an AFM cantilever. In Fig.3.1 
we have shown cartoons of three possible experimental setups. The two distinct procedures which 
can be followed to measure force- extension are: (a) Both the ends of the polymer are held via 
laser or magnetic tweezers, (b) One end of the polymer is attached to a substrate such that the 
position and orientation of this end is fixed while the other end is trapped via laser tweezer or 
magnetic tweezer or AFM cantilever. 

While optical tweezers allow free rotation of dielectric beads within the trap thereby allowing 
free orientations of the polymer end, magnetic tweezers fix the orientation of the ends and one 
can study the dependence of polymer properties on orientation of its ends by controlled change 
of the direction of external magnetic field . In this chapter, we call this fixing of orientation of 
an end of a polymer as grafting. By changing the trapping potential from stiff to soft trap one 
can go from Helmholtz to Gibbs ensemble[40]. Before we proceed, let us first elaborate on how 
to fix the ensemble of a mechanical measurement. In a simplest case we can assume that one 
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end of the polymer is trapped in a harmonic well, 

m^) - (3.1) 

with (0, 0, zq) being the position of potential- minimum. The polymer end will undergo continuous 
thermal motion. One can use a feedback circuit to shift zq to force back the fluctuating polymer 
end to its original position. This will ensure a Helmholtz ensemble. This can also be achieved by 
taking C ^ oo. On the other hand one can use a feedback circuit to fix the force — C(z — Zq) 
by varying C depending on the position z of the polymer end. This will ensure a Gibbs ensemble. 
This can also be achieved by taking a vanishingly soft (C — > 0) trap to infinitely large distance 
{zq oo) such that within the length scale of fluctuation the polymer end feels a constant slope 
of the parabolic potential. Surely, in experiments, using a feedback circuit is easier to implement 
a particular ensemble. However, the other procedure is mathematically well defined and one can 
seek recourse of it to show that the partition function of two ensembles are related by a Laplace 
transform [41]. This does not depend on the choice of Hamiltonian for the polymer. An exact 
relation between the two ensembles for worm like chain (WLC) model is shown in Sec. 3.1. 

From the above discussion on possible experiments, it is clear that there can be three possibili- 
ties of boundary conditions in terms of orientation. In an experiment the possibilities are, (a) free 
end: both the ends of a polymer can remain free to rotate[41, 42], (b) one end grafted: one end 
may be grafted and the other can take all possible orientations[43] and (c) both ends grafted: 
the orientations of both the ends are kept fixed. Thus, in experiments, we can have two possible 
ensembles and three possible boundary conditions. We investigate the probabilitiy distribution, 
free energy profile and force extension relation for each of these cases in this chapter. We shall 
see that the properties of a semiflexible polymer depend both on the choice of the ensemble and 
the boundary condition. 

WLC model is a simple coarse grained way to capture bending rigidity of an unstretchable 
polymer [21, 44] embedded in a thermal environment. Recent single molecule experiments in 
biological physics [18, 19, 26, 37] renewed interest in this old model of polymer physics. It 
was successfully employed by Bustamante, Marko, Siggia and Smith [25, 45] to model data 
from force- extension experiment [26] on double stranded DNA molecules. Mechanical properties 
of giant muscle protein titin [46, 47], polysaccharide dextrane [39, 47] and single molecule of 
xanthane [48] were also explained using WLC model. Due to the inextensibility constraint WLC 
model is hard to tract analytically except for in the two limits of flexible chain {t oo) and rigid 
rod {t 0), about which perturbative calculations have been done [27-29]. A key quantity that 
describes statistical property of such polymers is the end-to-end distance distribution. Numerical 
simulations for different values of t have been reported by Wilhelm and Frey [30], who have also 
obtained series expansions valid in the small t limit. Mean-field treatments by Thirumalai and 
his collaborators has also been reported [31, 49]. In an earlier study[42] we have investigated 
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the free energy profile of a semiflexible polymer whose ends were free to rotate in the constant 
extension ensemble and in the stiffness regime of 1 < i < 10. This has been predicted that 
a clear qualitative signature of semi- flexibility would be a non- monotonic force extension for 
stiffnesses around t — 4: m Helmholtz (constant end to end separation) ensemble. This comes 
from bimodality of probability distribution of end to end separation. This non- monotonicity must 
be absent in Gibbs (constant force) ensemble[42]. Multiple maxima in the probability distribution 
of end to end separation was due to a competition between entropy, that prefers a maximum near 
zero separation, and energy, that likes an extended polymer. A lot of theoretical works followed 
to reproduce and understand the probability density at all stiffnesses including the very interesting 
bimodality using analytic techniques [41, 50-52]. Recently, Frey and his collaborators studied the 
interesting multimodality in transverse fluctuations of a grafted polymer using simulations [43] 
and approximate theory [53, 54]. In a separate study Spakowitz and Wang used Greens function 
technique that takes into account the orientations of the polymer ends [55]. WLC model has 
been extended to study double stranded to single stranded DNA transition [56] and to incorporate 
twist degree of freedom [57-59]. 

The construction of this chapter is as follows. In Sec. 3.1 we present a path integral technique 
for exact calculation of WLC model for all the three boundary conditions and two ensembles. 
Then in Sec. 3. 2 we discuss the different discretized versions of WLC model and the Monte- Carlo 
(MC) simulation procedures followed in this work. In Sec. 4. 2 we present all the results obtained 
from theory and simulations. In this section we present statistici and mechanical properties for 
all the possible situations. Then we summarise our results and conclude with some discussions in 
Sec.5.7. 



3.1 Theory 

In WLC model the polymer is taken as a continuous curve denoted by a d- dimensional vector 
r(s) where s is a distance measured over the contour of the curve from any end of it. This curve 
has a bending rigidity and thus the Hamiltonian is given by 

where i{s) = df{s)/ds is the tangent vector and the polymer is inextendible i.e. i"^ = 1, P is the 
inverse temperature. The bending rigidity k is related to persistence length A via k = {d — l)X/2. 
Persistence length is a measure of the distance up to which the consecutive tangent vectors on 
the contour do not bend appreciably and is defined by (£(s).t(0)) = exp(— s/A). 

In this section we present a theoretical method to solve WLC model to any desired accuracy 
for both the Helmholtz and Gibbs ensemble and all three possible boundary orientations over the 
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entire range of stiffness parameter t. We first develop the method for a free polymer. Then we 
extend it to calculate properties of grafted [ one/both end(s) ] polymers. 

If the tangent vectors of two ends of a polymer are held fixed at ti and if, the probability 
distribution of end to end vector in constant extension ensemble can be written in path integral 
notation as 

P{r) ^-^ J^^^ T^iiis)] exp {-PH) x5^(^- tds^ (3.3) 

where PH is given by Eq.3.2 and I'[i(s)] denotes integration over all possible paths in tangent 
vector space from tangent at one end ii to tangent vector at the other end if. In d- dimensions 
r = (ri,r2, . . . ,rd). There is no eaxct analytic calculation of this distribution because of the 
difficulty presented by the inextensibility constraint introduced via the Dirac- delta function, 
though some mean field way of enforcing this constraint exist[31, 49]. Unfortunately that do 
not capture the very interesting triple maxima feature of the radial distribution at intermediate 
stiffness values as obtained in Ref.[42]. Recently, following our earlier work[42], J. Samuel et. 
al. [41] developed a path integral Greens function formulation to evaluate the distribution for a 
free polymer in 3D. We closely follow that method to generalize that to obtain results for various 
orientation constraints on polymer ends. 

The integrated (projected) probability distribution is given by, 

Px{x) ^ J dfP{r)5{ri- x). (3.4) 
We define the generating function of Px{x) via Laplace transform as, 

P{f) = J dxexp{Fx/kBT)Px{x) 

— J dxexp{fx/X)Px{x) (3.5) 
where / is the force in units of ksT/X i.e. / = FX/ksT applied along the x- axis. This gives, 

P(/) = M V[i{s)]e\ ' io'^Ha. j+AioW^j 
Jti 



f 



N / V[t{T')]e 



(3.6) 



The last step is obtained by replacing t' = s/\, t = L/ \ and using the identity k = {d— l)A/2. 
Note that, P{f), is the partition function in constant force ensemble where t behaves like an 
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inverse temperature such that the Gibb's free energy can be written as G{f) = —l/t lnP(/). 
Now considering r' as imaginary time and replacing r = —ir' we get, 

'■if 



P(f) ^ N I ' v[i{T)\ 



dT 



t, 

^ Af r V[i{r)]ei'^o'*Ldr] (3.7) 
Jii 

With the identification of L = ^ (%^)^ + /4 as the Lagrangian, P{f) [ = Z{f)/Z{0) ] in 
the above expression is the path integral representation for the propagator of a quantum particle, 
on the surface of a d- dimensional sphere, that takes a state \ti) to \tf). In Schrodinger picture 
this can be written as the inner product of a state and another state \if) evolved by imaginary 
time —it, 

Z{f) = {U\exp{-tH{-tt))\if) = {U\exp{-tH)\if). (3.8) 

Once P{f) — Z{f)/Z{0) is calculated, performing an inverse Laplace transform we can obtain 
the projected probability density Px{x). We now describe how to do that. Eq.3.5 can be written 
as, 

Pif) = J ^dxexp{fx/X)P,ix) 

= J dv^exp{tf V:c)LPx{x) 
"1 

dvxexp{tf Vx)px{vx) (3.9) 

-1 

where Vx = x/L and Px{vx) = LPx{x) is a scaling relation. Note that the Helmholtz free energy 
is given by J^xi^x) — —{^/'t)^'n.Px{vx)- Therefore, we can write the above expression as, 

exp[-tG(/)] = J dvx exp(t/ Vx) exp[-tJ^xivx)]. 

This is the relation between free energies of Helmholtz and Gibbs ensemble for a finite chain 
(finite t). In thermodynamic limit of t — > oo a steepest descent approximation of the above 
integral relation gives G{f) — Tx{vx) — fvx, the normal Legendre transformation. Let us use 
the identity, 

pAvx) = ^J e'^'^'^duj px{w)e-'^^dw (3.10) 
In it we can define Px{u) = ^\px{w) exp{—iuw)dw as the Fourier transform and 

Px{vx) = 7^ dupx{u) exp{iuvx) (3.11) 
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Figure 3.2: For a semiflexible polymer having its ends free to rotate Px{vx) ( = Py{vy) ) is plotted 
at stiffness parameter t = 2. The points are collected from Monte- Carlo simulation in freely rotating 
chain model (see Sec. 3. 2). The line is calculated from theory (see Sec. 3.1). The theory matches exactly 
with simulation. It clearly shows bistability via two maxima in integrated probability density at the two 
near complete extensions. 

as the inverse Fourier transform. Now, with —iu = tf we get P{f) = px{u = ift) and the 
inverse Fourier transform can be written as an inverse Laplace transform, 

Px{vx) = t— / dfPif) exp(-tM) (3.12) 
2vr2 y_ioo 

This gives the relation between the partition function P{f) in the constant force ensemble and 
the projected probability density of end to end separation Px{vx) along any given direction x in 
the constant extension ensemble. In numerical evaluation, the simplest way to obtain Pxiv^) is 
to replace / = —iu/t in the expression for Z(f) to obtain pxiu) and evaluate the inverse Fourier 
transform (Eq.3.11). 

Up to this point everything has been treated in d- dimensions. Experiments on single polymer 
can be performed in 3D as well as in 2D. In 3D polymers are left in a solution whereas one 
can float the polymer on a liquid film to measure its properties in 2D [16]. Moreover, polymers 
embedded in 2D are more interesting because of the following reason. It was shown that in 
Helmholtz ensemble in three dimensions[41], 

pK) = ^ (3.13) 

IttVx dVx 

where p{vx) is the scaled radial distribution function L'^P(r) = p{y) with v replaced by Vx- 
Since p{v) is a probability density p{v) > and therefore dpx/dvx < for t>a; > thus ruling 
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out multiple peaks in Px{vx) [40] and showing that Px{vx) will have single maximum at Vx = 
for all stiffness parameter t. No such simple relation exists between p{vx) and Px{vx) in two 
dimensions. Therefore, the two dimensional WLC polymer having its ends free to rotate may 
show more than one maximum in Px{vx) and therefore non-monotonicity. Indeed our calculation 
and simulation (see Sec.3.2) does show non-monotonicity (Fig. 3. 2). This is a curious difference 
between semiflexible polymers in 2D and 3D. Because of this and the fact that experiments in 
2D are possible, in this work we focus on the 2D WLC model. 

We have already given a general form of Z{f) (Eq.3.8) which depends on the dimensionality 
d of the embedding space of the polymer, d — 3 formulation of this theory for a polymer with 
its ends free to rotate has been carried out in detail in Ref.[41]. For d — 2, we can assume 
i — (cos ^, sin ^). Therefore, L — {1/4 9'^ + fcos9} and angular momentum pg — ^ — 9/2. 
Then the corresponding Hamiltonian is if = 9pe—L — pg—fcos9 and in planar polar coordinates, 
replacing pg ~^'§g we obtain the quantum hamiltonian operator, H — — ^ — fcos9. In this 
representation of t we can write, 

Zif) = {9,\expi-tH)\9f) 

= J2{W{n\^M-'tH)\n){n\9f) 

n,n' 

= J2€m<i>ni0f){n\exp{-tH)\n') (3.14) 

n,n' 

This propagator takes a definite tangent vector state \6i) at one end of the polymer to a definite 
final tangent vector state \9f) at the other end of it. 

AAA 

If external force is applied along x- direction as in Eq.3.5, H — Hq + Hj — —-^ — fcos9 
where Hq — — ^ is the Hamiltonian for a free rigid rotor in 2D and Hj — —fcos9 is the part 
of Hamiltonian introduced by an external field. Thus the total Hamiltonian H denotes a rigid 
rotor in presence of a constant external field. The eigenvalues of Hq, the hamiltonian for a 2D 
rigid rotor, are — and the complete set of ortho- normalized eigen- functions are given 

by 0„(^) = exp(m^)/-\/27r where n — 0, ±1, ±2, The orthonormality condition is {n\n') — 

(l/27r) J^'^ d9 ex.p[i{n - n')9] = Sn,n'- In this basis {n\Hi\n') = -{f /2){5n',n+i + Sn',n-i)- 
Therefore, {n\H\n') — T?bn\n — if /2){5n',n+i + (^n',n-i)- If the external force were applied in y- 
direction Hj = —fsm9 and {n\H\n') = n'^Sn',n — (//2i)(5n',n+i — <^n',n-i)- {n\ exp{—tH)\n') 
can be calculated by exponentiating the matrix {n\H\n'). Thus one can find Z{f). 

3.1.1 Free Polymer 

For a polymer which has both its ends free to rotate, integrating Eq.3.14 over all possible initial 
and final tangent vectors in rigid rotor basis we get 

Z{f) = 271 (0| exp(-t^)|0) (3.15) 



3.1 Theory 



25 



This means that Z{f) is given by the (0,0) element of the matrix {n\exp{—tH)\n') . In this 
case Z{0) = 2% (0| exp(-tifo)|0) = 2ti exp(-t 0^) = 27r and therefore P{f) = Z{f)/Z{Q) = 
(0| exp(— 1^/')|0). To evaluate the matrix element (0| exp(— ti^)|0) we exponentiate the matrix 
—t{n\H\n') and pick up the (0, 0)-th element. Thus, if the external force is applied in x- direction, 
remembering Px{u) = P{f = —iu/t) we calculate the inverse Fourier transform (Eq.3.11) to 
obtain Px{vx). Here it is useful to note that due to spherical symmetry of a polymer whose ends 
are free to rotate Px{vx) — Py{vy). 

3.1.2 One End Grafted 

This symmetry breaks down immediately if we hold one end of the polymer to a specific direction, 
namely along the x- direction i.e. 9i = 0. Then on Eq.(3.14) integrating over all possible 9f and 
leaving 6'^ = we obtain 

Z{f) = Y.{n\eM-tHm (3.16) 

n 

in the rigid- rotor basis. Note for this case Z{0) = ^„(n| exp(— tiJo)|0) = X^„exp(— tn^)5„^o = 
1 and therefore we have P{f) = Z{f). If the external force acts in x- direction, the Laplace 
transform of Z{f ), defined in the way described above, gives the projected probability distribu- 
tion in X- direction, Px{vx). On the other hand, if the external force acts in y- direction, the 
Laplace transform of Z{f) gives the projected probability distribution in y- direction Py{vy), the 
distribution of transverse fluctuation while one end of the polymer is grafted in x- direction. 

3.1.3 Both Ends Grafted 

Two ends of a polymer can be grafted in infinitely different ways. Let us fix the orientation of 
one end along x- direction {9i = 0) and the other end is orientated along any direction 9f, then 
Eq.(3.14) gives 

^(/) = ^ E ^^"'''(^1 eM-mn'). (3.17) 

n,n' 

As above, to obtain Px{vx) we use {n\Hi\n') — —{f /2){5n',n+i + Sn',n-i), whereas to obtain 

Py{vy) we use {n\Hi\n') — —{f /2i){Sn',n+i — (^n',n-i) to calculate P{f) and perform inverse 
Laplace transform. For this case of grafting 27rZ(0) = E„,n' = e''^^/"*'*" 
and therefore P(/) = X)n,n' e^"'^/(n|e-*^|n')/En e'"^^"*"']- 

Up to this point all the relations are exact. Since an infinite dimensional calculation of 
(n| cxp{—tH)\n') is not feasible, we calculate it numerically up to a dimension Nd, that controls 
the accuracy, limited only by computational power ^. We use = 10 which already gives 

^We use the MatrixExp function of Mathematica[60]. 
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numbers for probability distribution which differs within 1% from that obtained from = 20 
and gives a fair comparison with simulated data (see Sec.4.2). The inverse Laplace transform to 
obtain Px{vx) etc. from P{f) is also done numerically. 

3.2 Simulation 

In this section, we introduce two discretized models used to simulate semi- flexible polymers. 
Both of these are derived from the WLC model which is used for our analytical treatment in 
Sec.3.1. After introducing the discretized models we have listed the various boundary conditions 
used. We perform Monte- Carlo (MC) simulations to obtain probability distributions in Helmholtz 
ensemble. 

One discretized version of this model is the freely rotating chain (FRC) model[21]. In the FRC 
model one considers a polymer as a random walk of N steps each of length b — L/N with one 
step memory, such that, successive steps are constrained to be at an fixed angle 9 with A = 2b/9'^. 
The continuum WLC model is obtained in the limit 9, b ^ 0, N ^ oo keeping A and L finite. 
To simulate a polymer with ends free to rotate a large number of configurations are generated 
with first step taken in any random direction. Whereas if one choses the first step to be in some 
specific direction, this will simulate a polymer with one end grafted in that direction. 

A straight discretization of Eq.3.2 in 3D (2D) is an Id Heisenberg (classical xy) model: 

i=l i=l 

with a nearest neighbour coupling J = n/b between 'spins' ti. We have ignored a constant term 
in energy. The appropriate continuum limit is recovered for 6 — > 0, J — > oo with Jb — k finite. 
In this model grafting is simulated by fixing end spins on the ID chain. If an end is free then the 
end spin takes up any orientation that are allowed by the energy and entropy. In this model, by 
fixing the two end- spins, one can easily simulate a polymer with both its ends grafted in some 
fixed orientations. We follow the normal Metropolis algorithm to perform MC simulation in this 
model. 

We restrict ourselves to two dimensions. The numerics were checked via exact calculation 
of (r^) and (r^) which match within 0.5%. In the FRC model simulations we have used a 
chain length of = 10"^ and generated around 10*^ configurations. This simulation does not 
require equilibration run. Therefore all the 10*^ configurations were used for data collection. In 
xy model we have simulated N = 50 spins and equilibrated over 10^ MC steps. A further 10^ 
configurations were generated to collect data. We have averaged over 10'^ initial configurations, 
each of which were randomly chosen from nearly minimum energy configurations that conforms 
with the boundary conditions. Increasing the number of spins do not change the averaged data. 
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Figure 3.3: Compare this force- extension relation with that in Fig. 2. 4 of previous chapter. None of 
the curves including that at t = 3.33 show non- monotonic behaviour. This shows the force- extension 
behaviour in constant force ensemble. Forces are expressed in units of fc^T/A. 

3.3 Results 

Once all these theoretical and simulation tools are available, we apply them to bring out statistical 
and mechanical properties of a semiflexible polymer. We have three different boundary conditions 
depending on the orientations of polymer ends and two different ensembles. For each case we 
look at the various probability densities, ensemble dependence of force- extension etc. For the 
case of a polymer with both ends grafted we find that the properties depend on the relative 
orientation of the two ends. 

3.3.1 Free Polymer: 

Helmholtz ensemble: We employ the theory as described in Sec. 3.1 to calculate Pxivx) and Py{vy) 
for a polymer with both its ends free to rotate. We compare the probability distribution obtained at 
stiffness parameter t = 2 with that obtained from MC simulation (Sec. 3. 2) using the FRC model 
(see Fig. 3. 2). This shows exact agreement between theory and simulation. For a free polymer 
Px{vx) and Py{vy) are same due to the spherical symmetry. Note that, J^{vx) = —l/t\npx{vx) 
would give a non- monotonic force- extension via {fx) = (dJ-'/dvx). The force- extension obtained 
from the projected probability density Px{vx) will describe the experimental scenario in which the 
external potential traps the polymer end only in x- direction and constant in y. In general, if the 
external potential traps the polymer in dr dimensions {dr < d) then a dr dimensional projection 
( [d — dr) dimensional integration ) of the probability distribution of end to end vector gives 
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Figure 3.4: The simulation data for Px{vx) and Py{vy) from FRC model and XY- model simulations 
are compared with their theoretical estimates. Simulations and calculations were done at i = 2 for a 
polymer with one end grafted in x- direction. The data labelled LMF is taken from Ref.[43]. 

the free energy and decides the force- extension relation. This understanding is general and do 
not depend on the orientational boundary conditions at the polymer ends or the dimensionality 
of embedding space. This is important to keep in mind while analyzing experimental data. 

Gibbs ensemble: We have already mentioned that the non- monotonic nature of free energy, 
a strong quiitative signature of semi- flexibility, is observable only in Helmholtz ensemble not 
in Gibbs ensemble [42]. The partition function in Gibbs ensemble with force / applied in x- 
direction gives the Gibb's free energy G{f) = — ln[P(/)]. From this the averaged extension 
comes out to be {v^) = —{dG/df). For a polymer with its ends free to rotate, this force 
extension relations, that have been calculated from theory, at various t are shown in Fig. 3. 3. 
Notice that d{vx)/df = t[{vl) — (fx)^] > 0. Similar relation for response function does not exist 
in Helmholtz ensemble. Therefore, the force- extension in Gibbs ensemble has to be monotonic 
(Fig. 3. 3) in contrast to Helmholtz ensemble. Note that, for any non-zero stiffness, negative force 
is required to bring the end to end separation to zero. The amount of this force is larger for larger 
stiffness (smaller t). At large and positive force polymer goes to fully extended limit whichafter 
inextensibility constraint stops the polymer to extend any more. 

3.3.2 Grafted Polymer: One End 

Helmholtz ensemble: Next, we use our thery to plot Pxivx) and Py{vy) at t = 2 (Fig. 3. 4) for a 
semiflexible polymer with one end grafted in x- direction. In Fig. 3. 4 we have also plotted MC 
data from the FRC model and xy model simulation. The exact match validates our theory and 
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Figure 3.5: The integrated probability density Py{vy) is plotted at various stiffnesses t. At t = 4 
there is a single maximum at Vy = 0. Decreasing t we see at t = 2.8 emergence of two more peaks 
at nonzero Vy except for the one at = (See inset). At t = 2 central peak vanishes, the trimodal 
distribution becomes bimodal. The circles labeled LMF are data taken from Ref.[43] at f = 2 to show 
exact agreement with our theory. At t = 0.75 we see re-emergence of the central peak and tri-modality 
in Py{vy) (See inset, Os are from our MC simulation in the FRC model at t = 0.75, the lines are 
calculated from theory). 

both the simulation techniques. In Px{vx) the peak in near complete extension along positive x is 
due to the large bending energy and orientation of the other end towards this direction (also see 
p{vx,Vy) in Fig. 3. 13). We then explore the transverse fluctuation Py{vy) of this system, in detail, 
for different t (Fig. 3. 5). At large t{= 10), Py{vy) has single maximum at Vy = 0. At such low 
stiffnesses entropy takes over energy contributions. Number of possible configurations and thus 
entropy gains if end to end separation remains close to zero. This gives rise to this single central 
maximum. The emergence of multiple maxima at nonzero Vy in polymers, the multi- stability, 
with larger stiffness {t = 2.8) is due to the entropy- energy competition. The central peak is 
due to the entropy driven Gaussian behavior. The other two peaks emerge as entropy tries the 
polymer to bend around Vy = and energy restricts the amount of bending. Since bending in 
positive and negative y- directions are equally likely, the transverse fluctuation shows two more 
maxima at non- zero Vy. With further increase in stiffness {t = 2), the central Gaussian peak 
vanishes (see Fig. 3. 13) and therefore Py{vy) becomes bistable with two maxima (Fig. 3. 5). At 
even higher stiffness {t = 0.75) central peak reappears, this time due to higher bending energy. 
At t = 0.5 the distribution again becomes single peaked at f = as entropy almost completely 
loses out. This is the rigid rod limit. Notice that we have plotted MC data as obtained in Ref.[43] 
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Figure 3.6: The left panel shows the Helmholtz free energies J^{vx) and J^{vy) of a polymer with 
t = 2 and one end grafted in x- direction. The right panel shows the corresponding force- extensions 
in Helmholtz ensemble. Both {fx)- Vx and (fy)- Vy show regions of negative slope. Free energies are 
expressed in units of ksT and forces are expressed in units of ksT/X. 

for xy model simulation at t = 2. This agrees exactly with our theory. In the inset of Fig. 3. 5 we 
have blown up the multistability at t = 2.8 and t = 0.75. We have also plotted our FRC model 
simulation data at t = 0.75 for comparison. 

At this point it is instructive to look at the force extension behavior in Helmholtz ensemble, 
the ensemble in which Py{vy) and Px{vx) have been calculated above. In it the extension Vx [ or 
Vy ] is held constant and the corresponding average force in x- [ or y- ] direction is found from 
the relation (fx) = d!F{vx) / dvx ( or {fy) = dJ-'{vy) /dvy ). In Fig. 3. 6 we show the Helmholtz 
free energies J^{vx) = — (l/t)ln Px{vx) and J^{vy) = — (l/t)ln Py{vy) and the corresponding 
force extension curves in constant extension ensemble. Note that unlike the monotonicity beared 
by {vy)-fy curve (Fig. 3. 7) obtained in Gibbs ensemble the {fy)-Vy curve in Fig. 3. 6 clearly shows 
non-monotonicity, a signature of the Helmholtz ensemble. 

Gibbs ensemble: From our theory we can also explore the transverse response of a polymer 
which has one of its ends grafted to a substrate and a constant force is applied to the other 
end in a direction transverse to the grafting. Assume that the grafting direction is x and we 
apply a force fy in y- direction to study the response. A linear response theory was proposd 
earlier[53] to tackle this qustion. Our theory can predict the effect of externally applied force 
fy of arbitrary magnitude on the average positions {vx) and {vy). As the force is applied in 
^/-direction i.e. / = yfy, we have Hj = —fySinO. Because one end of the polymer is grafted 
in X- direction we use {n\Hj\n') = —{fy/2i){6n'.n+i — Sn',n~i) to evaluate Z{fy), whereas to 
calculate {vx) = —{dG/dfx) [ or, {vy) = —{dG/dfy) ] we introduce a small perturbing force 
5 fx ( or, 6fy ) in the Hamiltonian matrix to obtain the partial derivatives. Thus we obtain the 
corresponding force- extensions shown in Fig. 3. 7. As the grafted end is oriented in x- direction 
we expect in absence of any external force (f^) will be maximum and will keep on reducing due 
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Figure 3.7: Average displacements along x- direction {v^) and y- direction {vy) as a function 
of transverse force (transverse to grafting direction) in constant force ensemble. Lines denote our 
theoretical calculation while points denote the MC simulation data taken from Ref.[43]. Forces {Fy) are 
expressed in units of ksT/X, i.e. fy = FyX/kBT. 

to the bending of the other end generated by the external force imposed in y- direction. Thus 
{vx) is expected to be independent of the sign of fy. Similarly {vy) should follow the direction 
of external force and therefore is expected to carry the same sign as fy. Fig. 3. 7 verifies these 
expectations and shows the comparison of our theory with simulated data taken from Ref.[43]. 

3.3.3 Grafted Polymer: Both Ends 



Helmholtz ensemble: Let us first fix the orientations of the polymer at both its ends along x- axis 
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Figure 3.8: The simulation data for Px{vx) and Py{vy) from XY- model simulations of a WLC 
polymer are compared with their theoretical estimates. Simulations and calculations were done at t = 2 
for a WLC polymer with both its ends grafted in x- direction. 

and compare Px{vx) and Py{vy) obtained from our xy model simulation and our theory (Fig. 3. 8). 
This validates both our theory and simulation. Next, we go on to explore the properties of this 
system using our theory. We take a polymer with both its ends grafted, orientation of first end 
is fixed in x- dirction {9i = 0) and that in the other end {9f) is varied to study the change in 
transverse fluctuation Py{vy). To begin with let us find Py{vy) for different stiffnesses t with 
Of = (Fig. 3. 9). The height of central peak shows non- monotonicity; with increase in t from 
t = 1 the height of central peak first decreases up to t = 2 and then eventually again increases. 
The initial decrease in the height of the maximum is due to the fact that with increase in t the 
other end of the polymer (relative to the first end) starts to sweep larger angles about x- axis. 
With further increase in t entropic contributions that favor v = win over to increase the height 
of the maximum (see Fig. 3. 13). Note that, Py{vy) does not show multimodality as has been seen 
in the transverse fluctuation of a polymer with one end grafted. 

We then fix the stiffness at t = 4 and rotate the orientation of the other end and find out 
the transverse fluctuation Py{vy) (transverse to the orientation of first end, see Fig. 3. 9). At 
9 = Of = the fluctuation is unimodal with the maximum at Vy = 0. With increase in the 
orientation of the other end rotates from positive x- axis towards positive y- axis. Energetically the 
polymer gains the most if it bends on the perimeter of a circle. Therefore, energetically at any 0, 
the peak of py(vy) should be at f = ( 1. — cos6' )/0. Thus at = 0, 7r/4, tt/2, Sn/A, Tvr/S, vr the 
peak of py{vy) should be, respectively, at Vy = 0,0.37,0.64,0.72,0.69,0.64. Fig. 3. 9 shows that 
the peak positions almost follow these values up to 6* = tt/2, above which entropic contribution 
dominates. However, entropy always play a crucial role, e.g. at = n/A showing a double 
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Figure 3.9: The upper panel shows Py{vy) response of a polymer with both ends grafted along 
the same direction at various stiffness parameters t. They always remain single moded. In lower panel, 
Py{vy) plotted for various relative angle between the orientations of two ends at t = 4. The inset 
shows imergence of bistability at = 7r/4. 

peak around Vy = 0.37. At 6* = vr the two ends of the polymer are anti- parallel. Notice that, 
as = TT and 9 = —n are physically the same thing, Vy = ±0.64 are equally likely. Entropy 
would like the two ends to bend to Vy = 0. Competition between energy and entropy leads 
to almost a constant distribution up to \vy\ ~ 0.5. Notice that the bending energy couples to 
end orientations. For end orientations perpendicular to each- other, this coupling is strongest 
leading to a peak position that is closest to energetic expectation. The behavior of Py{vy) for 
—TV < 6* < is mirror symmetric about Vy = with respect to the behaviour of Py{vy) in the 
region < < vr. 

Gibbs ensemble: We then work in the constant force ensemble by applying a force / = yfy. Let 
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Figure 3.10: The left panel shows the transverse response and the right panel shows the longitudinal 
response of a polymer with both ends grafted along the same direction. 

us fix 9 = 0, i.e. both ends are oriented along x- axis. We find out the corresponding transverse 
and longitudinal response to this force in the similar manner as in the case of a polymer with 
one end grafted (see Fig. 3. 10). The force extension carries the same qualitative features as for 
a single end grafted polymer. For other end orientations {9 ^ 0) (v^) and (vy) does not remain 
at zero at fy = 0. Otherwise the nature of responses remain the same. To see the impact of 
the change in relative angle 9 between the orientations at the two ends we calculate the average 
extension of the polymer ends in x- and y- directions as we vary 9 in constant force ensemble 
at / = 0. Due to the bending energy, the {v^) is expected to by highest at = and lowest 
for 9 = ±7r. Similarly, (vy) is expected to be highest at = ± n/2 and lowest for 6* = and 
9 = ±7r. Fig. 3. 11 shows, while {vx) — 9 bears out the above expectations for all t, the expectation 
of getting highest magnitude for {vy) at 6* = ± 7r/2 is fulfilled only for lower stiffnesses [t > 5). 
At lower t the highest {vy) shifts towards higher relative angles 9 as bending energy takes over 
entropy. Note that, the smaller t gets, larger gets {vy), an impact of larger stiffness. 

3.3.4 End to End Vector Distribution 

We now employ MC simulations to study some other aspects of probability distribution in 
Helmholtz ensemble. We first examine the radial distribution function p(v) {v = r/L, p{v) = 
L^P{r)). It is clear from Fig. 3. 12 that grafting one end does not change the double maxima 
feature in p{v) at intermediate values of stiffness (4 < t < 2). Though grafting both the ends 
change the distribution function, the double maxima feature persists. Here we note that once 
one end of a polymer is grafted immediately the system loses its spherical symmetry, more so, 
since we restrict ourselves to semi-flexible regime. We have already seen that the projected prob- 
ability distribution Pxiv^-) and Py{vy) are very different depending on the orientations of polymer 
ends. However, for grafting of one end, this does not show up in radial distribution function p{v) 
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Figure 3.11: The upper panel shows the variation of {vx) as a function of 9 and the lower panel 
shows the variation of {vy) as a function of 6. {vx) and {vy) are calculated for stiffness parameters 

t = 1,2,3,4,5,10. 



(Fig. 3. 12). This is because, fixing one end only shifts the probability weight distributed over all 
possible angles at a given radial distance v towards the direction of orientation but on the same 
radial distance. 

The full statistics is encoded in the vector distribution function p{vx,Vy). To see the complete 
structure of it we next obtain p{vx,Vy) from MC simulations in FRC (for one or both ends of 
the polymer free to rotate) and xy- model (for polymer with both its ends fixed in some given 
orientation) and plot the total distribution p{vx,Vy) as a two dimensional density plot. We 
compare p{vx,Vy) of free polymer, polymer with one end grafted and polymer with both ends 
grafted (Fig. 3. 13). For definiteness we chose all the graftings to be in the x- direction. We 
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Figure 3.12: The radial distribution function p{v) at stiffness t = 4 is plotted for the three different 
boundary conditions - (a) both ends free, (b) one end oriented in x- direction and the other is free, (c) 
both ends oriented in x- direction. Radial distribution of first two cases are equal, whereas for the third 
case it is different. However, all three show double maxima. 

plot p{v3;, Vy) over a range of stiffnesses [t = 0.5, 2, 4, 10). The distribution has finite values for 
V < I and zero for v > 1. This is due to the inextensibility constraint in WLC model. In these 
density plots high probability is shown in red and low in blue (Fig. 3. 13). At small stiffness (large 
t) p{vx,Vy) shows a single entropic peak, at v = for free polymer and slightly shifted towards 
the direction of end- orientations in grafted polymers. This shifted entropic peak slowly moves 
to w = in the limit t oo. 

With increase in stiffness {t = 4) an energy dominated probability peak starts to appear near 
the full extension limit t> = 1 of the polymer. This peak forms a circular ring for free polymers. 
For a grafted polymer this peak is aligned in the direction of grafting. The probability distribution 
p{vx, Vy) at t = 4 clearly shows two regions of probability maxima, one near zero extension another 
near full extension, that gives rise to the double maxima in radial distribution function [42]. 

At larger stiffness (t = 2) the entropic maximum near the centre {v = 0) disappears. For 
the free polymer, one energy dominated maximum gets equally distributed over all angles. This 
way the system uses its spherical symmetry to gain in entropy. For grafted polymers, probability 
maximum near the full extension fans a finite solid angle around the direction of grafting. The 
distribution around the grafting direction is narrower for the polymer with both its ends grafted 
along the same direction. End orientation and bending energy couples to decide the probability 
distribution. This coupling is higher at smaller t [55]. This fact is more pronounced in plot of 
p{vx,Vy) at t = 0.5 (Fig. 3. 13). Spakowitz et. al. have investigated p{vx,Vy) using their Greens 
function calculations for a polymer with one of its ends grafted [55]. In this section we have used 
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Figure 3.13: Density plot of p{vx,Vy). Color code : red - high density, blue - low density. Left 
panel is for free polymer, middle panel is for polymer having one end grafted in x- direction and the 
right panel is for polymer having both ends grafted in x- direction. From top to bottom four panels 
denote stiffness values t = .5,2,4, 10. Note that the double maxima feature in p{vx,Vy) (one maximum 
near the centre and another near the rim) at t = 4 persists for all three boundary consitions. 

MC simulations to study p{vx,Vy) for all the possible boundary conditions. We have shown that 
multiple maxima in p{vx,Vy) persists near t = 4 for all these boundary conditions. 

As mentioned earlier, in Helmholtz ensemble the free energy is given by ^{vx, Vy) = — (1/t) \n[p{vx, Vy)]. 
In Fig. 3. 14 we plot this free energy profile on constant Vx = line at t = 4 and compare the three 
different boundary conditions. This is obtaind from Vx = cut of 2D probability density p{vx, Vy). 
This plot clearly shows that triple minima in free energy [42] prevails even after grafting one or 
both ends of a semi- flexible polymer. This means that in constant extension ensemble if one 
fixes the position of one end and traps the other end in 2D (both x- and y- directions) to extend 
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Figure 3.14: At i = 4 free energy profile J^{vy) in y- direction at = corresponding to the 
probability distributions shown in Fig. 3. 13 are plotted. This clearly shows that the tripple minima feature 
in free energy for a polymer with both ends free persists even after grafting one or both ends of the 
polymer. 



the end to end separation of a semiflexible polymer from the first free energy minimum at Vy = 
to full extension of Vy along Vx = line, for small extensions an average force will try to pull back 
the separation to zero corresponding to the central minimum of the free energy. However, once 
the separation is taken beyond the maxima in the free energy near Vy = ±0.5 the average force 
will push the separation away from the center towards the non- zero Vy minima in free energy. 
Further extension will cause a huge force towards this second minimum in the free energy. This 
is due to the inextensibility constraint in WLC model. Thus, in force- extension experiments on 
a polymer in constant extension ensemble, this multi- stability (non- monotonic force extension) 
at intermediate stiffness values will be measurable for all kind of boundary conditions[42]. The 
coupling of end orientation and bending energy have raised the free energy of a polymer with 
both its ends grafted, with respect to that of a free polymer and a polymer with one of its grafted 
(Fig.3.14). 

We again emphasize that whether the end to end vector distribution or some projected prob- 
ability distribution will capture the results of a force- extension measurement will depend on the 
kind of trapping potential (what are the directions in which it traps a polymer end). Apart from 
this, as we have shown, the orientational boundary conditions at the ends of a polymer and the 
ensemble of experiment will affect the force- extension behavior non- trivially. 
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In this chapter, we found the exact solution of WLC model of semi- flexible polymers for all 
possible end orientations and in both Gibbs and Helmholtz ensembles. In previous chapter, we 
had emphasized on the ensemble dependence of mechanical properties [42]. In this chapter we 
have shown that polymer properties also crucially depend on end orientations. Conditions of fixed 
end orientations are biologically important in several processes. Magnetic tweezer experiments 
can be used in laboratory to study orientation dependent properties of semi- flexible polymers. 
We have presented an exact theory to obtain end to end distribution function of WLC polymers by 
calculating the propagator for a quantum particle moving on a unit sphere. For a finite chain, the 
free energies in two ensembles were shown to be related via a Laplace transform relation which 
in thermodynamic limit goes over to the typical Legendre transform relation. Two discretized 
versions of WLC model have been used in MC simulations. In Helmholtz ensemble, for stiffness 
parameters around t = 4, we found multimodality (multiple minima in free energy) for all kinds 
of boundary orientations. Thus we have generalized the finding of triple minima in free energy 
that we obtained for free polymers in the previous chapter [42]. The radial distribution function 
of a free polymer and a polymer with a fixed orientation in one end show the same behavior. 
This we trace back to the properties of full distribution of end to end vector. We have calculated 
various projected probability distributions. The transverse fluctuation of a polymer with a fixed 
orientation at one end had been studied using MC simulation [43] and an effective medium theory 
[54] which approximately captures the qualitative features of the simulations. We have calculated 
this fluctuation from our theory that shows exact numerical agreement with our simulation and 
Lattanzi et. al.'s [43] simulation results. For a polymer with orientations at both its ends fixed, 
the end to end distribution shows dependence on relative angle 9 between the two ends. The end 
orientations couple to the bending energy to take an unimodal distribution through bimodality 
with change in 9. We have seen the general presence of non- monotoinic force extension in 
Helmholtz ensemble. We have also studied the force- extension behavior in Gibbs ensemble. 
These never show non- monotonicity. It is important to note that, multiple maxima in end to 
end distribution means a non- monotonic force extension in Helmholtz ensemble, not in Gibbs 
ensemble. As we have shown, the response function in Gibbs ensemble is always positive, thus non- 
monotonicity in this ensemble is impossible. We have outlined the general theoretical framework 
in d- dimensions. However, in this study we have looked at the properties of a 2D polymer. At the 
onset we have shown that an important point of polymer statistics, multimodality, is dependent 
on the dimensionality of embediing space. In three dimensional free polymer, multimodality in 
projected probability distribution is impossible, however presence of this is a reality in 2D. Similar 
studies in 3D remains to be an interesting direction forward. The presence of multimodality in 
intermediate range of stiffness (near t = 4) gives rise to a possibility of multi-stability in the 
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dynamics of semi- flexible polymers. The simplest question remains, for a free polymer how does 
the presence of a triple minima in free energy affects its dynamics. This might be important in 
the context of understanding the time- scale of a messanger- RNA finding out an active site on 
a DNA chain. 



4 Laser Induced Freezing and Re-entrant Melting 



Reviewing Melquiades' notes, serene now, without the exaltation of novelty, in 
prolonged and patient sessions they tried to separate Ursula's gold from the 
debris that was stuck to the bottom of the pot. - G. G. Mdrquez 



In the previous two chapters we have studied the properties of a finite semiflexible polymer in 
different ensembles and under different conditions of its boundary orientations. These conditions 
can be implemented on the polymer by imposing various kinds of external trapping and constraints. 
The finite size and bending rigidity of such molecules coupled with the boundary orientations 
rendered the different statistical and mechanical properties. In this chapter we consider another 
kind of trapping. A one dimensional modulating external trapping potential that couples with the 
local density is applied on a two dimensional system of particles interacting via various spherically 
symmetric short range repulsions. This trapping potential constrains the system in one direction. 
We find out the phase behavior of this system with change in the strength of the trapping 
potential. It is known that in two dimensions melting occurs by the unbinding of dislocation pairs 
[61-63]. How does this scenario change if there are periodic constraining potentials which may 
reduce the dimensionality of the system from two to one? 

Examples of phase transitions mediated by the unbinding of defect pairs abound in two di- 
mensions. The quasi- long- ranged order to disorder transition in the XY and planar rotor 
models[61, 64-68], the melting transition of two dimensional solid[61-63, 69-71], the super- 
conductor to normal phase transition in two dimensional Josephson junction arrays[72], the 
commensurate- incommensurate transition of the striped phase of smectic liquid crystals on 
anisotropic substrates[73], and the more recent discovery of a defect mediated re-entrant freezing 
transition in two dimensional colloids in an external periodic potential [74, 75] are all understood 
within such defect unbinding theories. While the very first defect mediated transition theory 
for the phase transition in the XY-modei by Kosterlitz and Thouless (KT)[61] enjoyed almost 
immediate acceptance and was verified in simulations[64-66, 68] as well as experiments[76, 77], 
defect mediated theories of two dimensional melting took a long time to gain general acceptance 
in the community[78]. There were several valid reasons for this reticence however. 

Firstly, as was recognized even in the earliest papers[62, 63] on this subject, the dislocation 
unbinding transition, which represents an instability of the solid phase, may always be pre-empted 
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by a first order[79, 80] transition from a metastable solid to a stable liquid. Whether such a first 
order melting transition actually occurs or not depends on the temperature of instability Tkt', so 
that if the transition temperature Tc < Tkt the unbinding of dislocations does not occur. Clearly, 
neither this condition nor its converse can hold for all 2d systems in general. This is because 
Tkt is a non-universal number which depends on the "distance" in coupling parameter space 
between the bare and the fixed point Hamiltonian and hence on the details of the interaction. 
Secondly, the renormalization group flow equations derived in all defect mediated theories to date 
are perturbative expansions in the defect density (fugacity) in the ordered phase. How fast does 
this perturbation series converge? The answer lies again in the position of the bare Hamiltonian 
in the coupling parameter space. For the planar rotor model[61, 68], past calculations show 
that next to leading order terms in the flow equations are essential to reproduce the value of the 
transition temperature obtained in simulations[68]. Thirdly, defect mediated transitions predict an 
essential singularity[61] of the correlation length at the transition temperature. This means that 
effects of finite size[81] would be substantial and may thoroughly mask the true thermodynamic 
result. A rapid increase of the correlation length also implies that the relaxation time diverges 
as the transition temperature is approached - critical slowing down. For the two dimensional 
solid, this last effect is particularly crucial since, even far from the transition, the motion of 
defects is mainly thermally assisted and diffusional and therefore slow. The equilibration of 
defect configurations[82] is therefore often an issue even in solids of macroscopic dimensions. 

On the other hand, over the last few years it has been possible to test quantitatively some of 
the non-universal predictions of defect mediated theories of phase transitions using simulations 
of restricted systems[68, 71, 83]. A simulation of a system without defects is used to obtain 
the values for the bare coupling constants which are then taken as inputs to the renormalization 
group equations for the appropriate defect unbinding theory to obtain quantities like the transition 
temperature. Needless to say, the simulated system does not undergo a phase transition and 
therefore problems typically related to diverging correlation lengths and times do not occur. 
Numerical agreement of the result of this calculation with that of unrestricted simulations or 
experiments is proof of the validity of the RG flow equations[61-63, 71]. This idea has been 
repeatedly applied in the past to analyze defect mediated phase transitions in the planar rotor 
model[68], two dimensional melting of hard disks[71] and the re-entrant freezing of hard disks in 
an external periodic potential [83, 84]. The last system is particularly interesting in view of its 
close relation with experiments on laser induced re-entrant freezing transition in charge stabilized 
colloids [74, 75] and this constitutes the subject of the present chapter as well. 

In this chapter we show in detail how restricted simulations of systems of particles interact- 
ing among themselves via a variety of interactions and with a commensurate external periodic 
potential can be used to obtain phase diagrams showing the re-entrant freezing transition. The 
results obtained are compared to earlier unrestricted simulations for the same systems. Briefly 
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our results are as follows. Firstly, we observe that, as in an earlier study of the planar rotor 
model[68], next to leading order corrections to the renormalization flow equations are essential 
to reproduce even the gross features of the phase diagram. Specifically, the re-entrant portion 
of the phase diagram can be reproduced only if such correction terms are taken into account. 
Secondly while we find almost complete agreement with earlier results for the hard disk system 
which has been studied most extensively, our phase diagram for the other forms of interaction is 
shifted with respect to the results available in the literature. This may mean either of two things 
— inadequacy of the RG theory used by us or finite size effects in the earlier results. Lastly, as 
a by product of our calculations, we have obtained the core energy for defects (dislocations) in 
these systems and studied its dependence on thermodynamic and potential parameters. 

The problem of re-entrant freezing transition of a system of interacting colloidal particles in a 
periodic potential has an interesting history involving experiments[74, 75], simulations[85-92] and 
theory[93-95]. In last couple of decades soft systems like colloids have been studied extensively[96] 
both for their own sake and as typical toy models to study various important condensed matter 
questions like structural and phase transitions through experiments that allow real space imaging. 
Charged colloids confined within two glass plates form a model 2-d system as the electrostatic 
force from the plates almost completely suppresses the fluctuations of colloids perpendicular to the 
plates, practically confining them to a 2-d plane. In their pioneering experiment Chowdhury[74] 
et. al. imposed a simple static background potential which is periodic in one direction and 
constant in the other (except for an overall Gaussian profile of intensity- variation) by interfering 
two laser beams. This potential immediately induces a density modulation in the colloidal system. 
The potential minima are spaced to overlap with the close packed lines of the ideal lattice of 
the colloidal system at a given density. With increase in potential strength such a colloidal 
liquid has been observed to solidify. This is known as laser induced freezing (LIF). In a recent 
experiment[75] it has been shown that with further increase in potential strength, surprisingly, 
the solid phase re-melts into a modulated liquid. This phenomenon is known as re-entrant 
laser induced freezing (RLIF). Qualitatively, starting from a liquid phase, the external periodic 
potential immediately induces a density modulation, reducing fluctuations which eventually leads 
to solidification. Further increase in the amplitude of the potential reduces the system to a 
collection of decoupled 1-d strips. The dimensional reduction now /ncreases fluctuations remelting 
the system. 

The early mean field theories, namely. Landau theory[74] and density functional theory[93] pre- 
dicted a change from a first order to continuous transition with increase in potential strength and 
failed to describe the re-entrant behavior, a conclusion seemingly confirmed by early experiments[74] 
and some early simulations[85]. Overall, the results from early simulations remained inconclusive 
however, while one of them[85] claimed to have found a tri-critical point at intermediate laser in- 
tensities and re-entrance, later studies refuted these results [86-88]. All of these studies used the 
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change in order parameter and the maximum in the specific heat to identify the phase transition 
points. While some later studies[86-88] found RLIF for hard disks they reported laser induced 
freezing and absence of any re-entrant melting for the DLVO potential[88] in direct contradiction 
to experiments[75]. 

Following the defect mediated disordering approach of Kosterlitz and Thouless [61](KT), Prey, 
Nelson and Radzihovsky[94](FNR) proposed a detailed theory for the re-entrant transition based 
on the unbinding of dislocations with Burger's vector parallel to the line of potential minima. 
This theory predicted RLIF and no tricritical point. The results of this work were in qualitative 
agreement with experiments[75] and provided a framework for understanding RLIF in general. 
More accurate simulation studies on systems of hard disks[89], soft disks[91, 92], DLVO[90] etc. 
confirmed the re-entrant freezing-melting transition in agreement with experiments[75] and FNR 
theory[94, 95]. In these studies the phase transition point was found from the crossing of Binder- 
cumulants[5, 97] of order parameters corresponding to translational and bond- orientational order, 
calculated for various sub- system sizes. A systematic finite size scaling analysis[89] of simulation 
results for the 2-d hard disk system in a 1-d modulating potential showed, in fact, several universal 
features consistent with the predictions of FNR theory. It was shown in these studies that the 
resultant phase diagram remains system size dependent and the cross- over to the zero field 
KTHNY melting[62, 63] plays a crucial role in understanding the results for small values of 
the external potential. While the data collapse and critical exponents were consistent with KT 
theory for stronger potentials, for weaker potentials they match better with critical scaling[89]. A 
common problem with all the simulation studies might be equilibration with respect to dislocation 
movements along climb (or even glide) directions. Also, non universal predictions, namely the 
phase diagram are difficult to compare because the FNR approach (like KT theory) is expressed 
in terms of the appropriate elastic moduli which are notoriously time-consuming to compute near 
a continuous phase transition. Diverging correlation lengths and times near the phase transition 
further complicate an accurate evaluation of the non universal predictions of the theory. 

We calculate the phase diagrams of three different 2-d systems with a 1-d modulating potential 
(see Fig. 6.1) following the technique of restricted Monte Carlo simulations[68, 71, 83], to be 
discussed below. For the laser induced transition we use this method to generate whole phase 
diagrams. We reject Monte Carlo moves which tend to distort an unit cell in a way which changes 
the local connectivity[71]. The percentage of moves thus rejected is a measure of the dislocation 
fugacity[71]. This, together with the elastic constants of the dislocation free lattice obtained 
separately, are our inputs (bare values) to the renormalization flow equations[94, 95] to compute 
the melting points and hence the phase diagram. Our results (Fig. 4.13,4.17,4.16) clearly show a 
modulated liquid (ML) locked floating solid (LFS) ML re-entrant transition with increase 
in the amplitude (Vq) of the potential. In general, we find, the predictions of FNR theory to be 
valid. 
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Figure 4.1: This cartoon shows a typical 2-d system under consideration, d is the length scale 
over which repulsive two body potentials are operative. The dashed lines indicate minima of external 
modulating potential l3V{y) = — /?Vb cos(27ry/(io). ax = ao is the lattice parameter fixed by the density 
p and Qy indicate the average separation between two layers along y-direction perpendicular to a set 
of close-packed planes. For a perfect triangular lattice ay = \/3ao/2. The modulating potential is 
commensurate with the lattice such that do = ay. 

Lastly, we must mention that our technique, as summarized above and used in this as well 
as earlier[83] work, corresponds closely with early studies of the melting of the electron solid by 
Morf[69, 98]. In the latter case, the dislocation fugacity, which is one of the important inputs 
to the KTHNY flow equations was obtained by a careful and direct calculation of the dislocation 
core energy at T = 0. Our approach is somewhat cruder but gives us numbers for nonzero T 
which automatically contain the effects of phonon fluctuations. 

In section 4.1 we first briefly discuss the FNR theory and then go on to show in detail the 
restricted simulation scheme used by us to obtain the various quantities required to calculate the 
phase diagram. In section 4.2 we give the simulation results. We describe, in detail, the various 
quantities leading to the phase diagram for one of the systems, viz. the hard disks[71, 99]. Then 
we present the phase diagrams for the other two systems we study. We compare our results with 
earlier simulations. Lastly, in section 5.7 we summarize our main results and conclude. 

4.1 Method 

A cartoon corresponding to the systems considered for our study is given in Fig. 6.1. The elastic 
free energy of the solid is given in terms of the spatial derivatives of the displacement field 
u{r) = r — fo with tq being the lattice vectors of the undistorted reference triangular lattice. 
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For a solid in presence of a modulating potential /3V{y) (Fig. 6.1) the displacement mode Uy 
becomes massive, leaving massless Ux modes. After integrating out the Uy modes the free energy 
of the LFS may be expressed in terms of gradients of and elastic moduli[94, 95], namely, the 
Young's modulus K{l3Vo,p) and shear modulus n{(3Vo,p), 



Hei — / dxdy 



1 / 1 / du. ^ ^ 



(4.1) 



Similar arguments[94, 95] show that among the three sets of low energy dislocations available in 
the 2-d triangular lattice, only those (type I) with Burger's vector parallel to the line of potential 
minima survive at large PVq. Dislocations with Burger's vector pointing along the other two 
possible close-packed directions (type II) in the 2-d triangular lattice have larger energies because 
the surrounding atoms are forced to ride the crests of the periodic potential[94, 95]. Within this 
set of assumptions, the system therefore shares the same symmetries as the XY model. Indeed, 
a simple rescaling of x — > ^y]Ix and y y/Ky leads this free energy to the free energy of the 
XY-model with spin-wave stiffness K^y — ^/Kjlal/An'^ and spin angle 9 — 27rUx/ao: 



Hei = / dxdy 



l^xyiver 



This immediately leads to the identification of a vortex in XY model (^ d9 — 2tt) with a dislocation 
of Burger's vector b — iao (^ du^ — ao, i — unit vector along x- direction) parallel to the potential 
minima i.e. the dislocation of type I. The corresponding theory for phase transitions can then be 
recast as a KT theory[61] and is described in the framework of a two parameter renormalization 
flow for the spin-wave stiffness (3Kxy{l) and the fugacity of type I dislocations y'{l), where I is a 
measure of length scale as / = ln(r/ao), r being the size of the system. The flow equations are 
expressed in terms of x' = {irKxy — 2) where K^y = PK^y and y' = Att exp{—l3Ec) where E^, is 
the core energy of type I dislocations which is obtained from the dislocation probability[71, 98]. 
Keeping up to next to leading order terms in y' the renormalization group flow equations[68, 100] 
are. 



dx' 

m 



Flows in / generated by the above equations starting from initial or "bare" values of x' and y' 
fall in two categories. If, as / ^ oo, y' diverges, the thermodynamic phase is disordered (i.e. 
ML), while on the other hand if y' vanishes, it is an ordered phase (a LFS)[94, 95]. The two 
kinds of flows are demarcated by the separatrix which marks the phase transition point. For 
the linearized equations, that keeps i only the leading order terms in y', the separatrix is simply 
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the straight line y' = x' , whereas for the full non-linear equations one needs to calculate this 
numerically[68, 71, 100]. 

The bare numbers for x' and y' are relatively insensitive to system size since our Monte Carlo 
simulation does not involve a diverging correlation length associated with a phase transition. This 
is achieved as follows[68, 71]. We monitor individual random moves of the particles in a system 
and look for distortions of the neighboring unit cells. If in any of these unit cells the length 
of a next nearest neighbor bond becomes smaller than the nearest neighbor bond, the move is 
rejected. All such moves generate disclination quartets and are shown in the Fig. 4.2. Notice 
that each of these moves break a nearest neighbour bond to build a new next nearest neighbour 
bond, in the process generating two 7-5 disclination pairs. These are the moves rejected in 
the restricted simulation scheme we follow. The probabilities of such bond breaking moves are 
however computed by keeping track of the number of such rejected moves. One has to keep track 
of all the three possible distortions of the unit rhombus with measured probabilities P^j, i = 1,2,3 
(see Fig. 4.2), 

^ number of rejected bond breaking moves of type i 
Total number of MC moves 
Each of these distortions involves four 7 — 5 disclinations i.e. two possible dislocation- antidislo- 
cation pairs which, we assume, occur mutually exclusively in a way that we explain shortly. For 
a free {Vq — 0) two dimensional system the dislocation core energy El can be found through the 
relation [98] 

n = exp(-/3 2£;*)Z(X) (4.3) 

where 11 = Yl^=iPmi ^he total number density of dislocation pairs per particle and Z{K) 
is the "internal partition function" incorporating all three types of degenerate orientations of 
dislocations, 

where Jq is a modified Bessel function, K = fiKa^ is a dimensionless Young's modulus renor- 
malized over phonon modes, being the lattice parameter and r„ij„ is the separation between 
dislocation-antidislocation above which one counts the pairs. The above expression for Z{K) 
and Eq.(4.3) have been used previously in simulations[71, 98] of phase transitions of 2-d systems 
in absence of any external potential to find the dislocation core energy E^. 

We now show how the probabilities for generating pairs of specific types of dislocations Pdi for 
Vo 7^ are related to P^i- Consider Fig. 4. 2 where each of the three varieties of bond breaking 
moves are depicted. It is clear that each given distortion can occur due to the presence of 
two possible dislocation- antidislocation pairs acting independently. For example a distortion of 
type- A can take place either due to dislocation dipoles with Burgers vectors making an angle 
of 60" with the horizontal or an angle of 120°. Both of these dislocation dipoles are of type II. 
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Figure 4.2: This diagram depicts all the possible dislocation generating moves that we reject. 
Starting from the triangular lattice shown in the centre (the dotted lines show the potential minima), 
in all, there can be three types of dislocation- pair generating moves shown as A, B & C. The numbers 
7 and 5 denote the positions of two types of disclinations having seven nearest neighbours and five 
nearest neighbours respectively. Only those bonds, which are necessary to show distortions due to the 
generation of disclination quartets, have been drawn. The rhombi near each of the distorted lattice 
denote the unit cells and open arrows from 7^5 show the direction of dislocation generating moves. 
The probabilities of these moves are Pmi(A), Pm2(B) and Pm3(C). Corresponding Burger's vectors 
(filled arrows) are bisectors pointing towards a direction rotated counter-clockwise starting from 7^5 
directions and are parallel to one of the lattice planes. Notice, the separation between Burger's vectors 
of a pair along the glide direction (parallel to the Burger's vectors) is a single lattice separation (oq) 
and within this construction it is impossible to draw Burger's loop that can generate non-zero Burger's 
vector. Depending on which of the two possible disclination pairs separate out any one dislocation- 
antidislocation pair will be formed. 

If this bond breaking move were to be accepted, then at a subsequent time step the individual 
dislocations making up any one of the two possible pairs could separate, the two possible events 
being mutually exclusive. This allows us to write down the following relations among the various 
probabilities. 

Pml — Pd2 + Pd3: Pm2 — Pd2 + Pdl: Pm3 — Pd3 + Pdl- 

Solving for P^jS and remembering that P^2 = Pds by symmetry we get Pai = l{Pr,i2 + Pm3-Pmi) 
and Pd2 = -Pmi/2. The above expressions are motivated and illustrated in Fig. 4. 2. Once the 
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probability of dislocation pairs are obtained in this fashion, we may proceed to calculate the 
dislocation core energy Ec and the dislocation fugacity y'. 

An argument following the lines of Fisher et. a/.[98] shows that the dislocation probability 
(number density of dislocation pair per particle) for our system, 

Pdi = exp(-/3 2E,)Z{K,y) (4.4) 

where 2Ec is the core energy and Z^K^y) is the internal partition function of dislocation pair of 
type I (single orientation). 



f (fr 
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(4.5) 



nK^y — 1 

with K^y — PKxy and Ac — y/3al/2 being the area of an unit cell in the undistorted lattice. 
We choose r^j„ = 2ao. At this point this choice is arbitrary. We give the detailed reasoning 
for this choice at the end of the discussions on hard disks in the section 4.2. Eq.4.4 and Eq.4.5 
straightaway yield the required core energy Ec. The corresponding fugacity contribution to RG 
flow equations (Eq.4.2) is given via 



y'^An^Pdi/Z{K,y) (4.6) 

In the above, the following assumption is, however, implicit. Once a nearest neighbor bond 
breaks and a potential dislocation pair is formed, they separate with probability one ^.This as- 
sumption goes into the identity Eq.4.4 as well as in Eq.4.3[71]. Taking the rejection ratios due 
to bond- breaking as the dislocation probabilities, as well, require this assumption ^. 

The same restricted Monte Carlo simulation can be used to find out the stress tensor, and the 
elastic moduli from the stress-strain curves. The dimensionless stress tensor for a free (Vq = 0) 
system is given by [102], 

NS,. I (4.7) 




where i, j are particle indices and A, u denote directions x, y; 0(r'^) is the two- body interaction, 
S/d"^ is the dimensionless area of the simulation box ^. 

"^This assumption is similar in spirit to assuming that a particle which reaches the saddle point in the 
Kramers barrier crossing problem would automatically cross the barrier [101]. 

^Note that the calculation of the bare fugacity from the dislocation probability is, we believe, more accurate 
that the procedure used in [83]. 

^In the presence of an external ID modulating potential periodic in the y- direction the stress has contri- 
bution from another virial- like additive term, — ^g- V^fy)' where is the y-component of position 
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4.2 Results and Discussion 

In this section we present the results from our simulations for three different 2-d systems, namely 
hard disks, soft disks and a system of colloidal particles interacting via the DLVO (Derjaguin- 
Landau-Verwey-Overbeek)[103, 104] potentials. We discuss, first, our calculation for a two 
dimensional system of hard disks, in detail. The procedure followed in other systems is almost 
identical. 

Hard disks: The bulk system of hard disks where particles i and j, in 2-d, interact via the 
potential 0(r*^) = for r*^ > d and 0(r'-') = oo for r*^ < d, where d is the hard disk diameter 
and r*^ = jr-' — r*| the relative separation of the particles, is known to melt[71, 99, 105-107] 
from a high density triangular lattice to an isotropic liquid with a narrow intervening hexatic 
phase[62, 63, 71, 99]. The hard disk free energy is entirely entropic in origin and the only 
thermodynamically relevant variable is the number density p = N/V or the packing fraction 
77 = (7r/4)pd^. Simulations[99], experimental [96] and theoretical[108] studies of hard disks show 
that for 7] > .715 the system exists as a triangular lattice which transforms to a liquid below 
7] = .706. The small intervening region contains a hexatic phase predicted by the KTHNY 
theory[62, 63] of 2-d melting. Apart from being easily accessible to theoretical treatment[109], 
experimental systems with nearly "hard" interactions viz. sterically stabilized colloids[96] are 
available. 

In presence of a periodic external potential, the only other energy scale present in the system 
is the relative potential ^strength PVq. If the modulating potential is commensurate with the 
spacing between close- packed lines, the elastic free energy of this system in it's solid phase 
follows Eq.4.1 and the corresponding renormalization flow equations are given by Eq.4.2. 

We obtain the bare y' and x' from Monte Carlo simulations of 43 x 50 = 2150 hard disks and 
use them as initial values for the numerical solution of Eqs. (4.2). The Monte Carlo simulations 
for hard disks is done in the usual[110] way viz. we perform individual random moves of hard disks 
after checking for overlaps with neighbours. When a move is about to be accepted, however, 
we look for the possibility of bond breaking as described in the previous section (Fig. 4. 2). We 
reject any such move and the rejection ratios for specific types of bond breaking moves give us 
the dislocation probabilities of type I and II, separately (Fig. 5. 12). 

vector of particle A. This contribution comes from the part of the free energy that involves higher energy 
(massive) excitations. For the elastic free energy which is lowest order in the displacement gradient 
(Eq.4.1) this part docs not contribute towards the elastic constants, as the x- and y- component of gra- 
dient remain uncoupled. This extra term in stress remains a constant background without disturbing 
the elastic constants connected to the Young and shear modulus that corresponds to distortions of the 
system in the low energy directions. We therefore neglect this background in calculating stresses where 
from we obtain the elastic moduli. 
^This interaction in colloids is due to polarization of the dielectric colloidal particles by the electric field 
of the laser. Though experiments of Refs.[74, 75] use charged colloids, the interaction of hard sphere 
colloids with lasers is similar. 
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Figure 4.3: Number density of dislocation pairs of type I and II per particle as a function of the 
amplitude of the laser potential fSVo. In this plot the O symbols correspond to P^i. the probability 
for type I dislocations and the X symbols to Pd2 the probability for type II dislocations obtained 
from the Pmi (see text and Fig. 4. 2) for various rj values, arrow denoting the direction of increasing 
r]{= .69, .696, .7029, .71). The Pji for i = 1,2 are multiplied by 10^. These probabilities are plotted 
against the potential strength PVq. Note that for /5Vo > 1, the probability for type I dislocations is 
larger than that of type II. The dots and solid lines are only guides to eye. 

From Fig. 5. 12 it is clear that the probability of type II dislocations i.e. Pd2 drops down to 
zero for all packing fractions at higher potential strengths I3Vq. The external potential suppresses 
formation of this kind of dislocations. For small (3Vq on the other hand, the probabilities of type I 
and type II dislocations are roughly the same. This should be a cause of concern since we neglect 
the contribution of type II dislocations for all (5Vq. We comment on this issue later in this section. 
Using Eq.4.6 and Eq.4.5 along with the identity r^i„ = 2ao gives us the initial value to be 
used in renormalization flow Eq.4.2. 

Before we move on, we comment on the magnitude of the errors for Pmi and hence for 
^Q. There are two main sources of errors for these quantities. This may arise from (a)finite 
simulation times and (b)small size of the system. In order to check for this, we have plotted the 
accumulated values for the probability Pai as a function of Monte Carlo step for 2150 and 21488 
hard disks (Fig. 4. 4). It is clear that our estimates for the probabilities are virtually error free! 
This demonstrates clearly the usefulness of our restricted Monte Carlo scheme. 

To obtain K^y we need to calculate the Young's modulus K and shear modulus In order to 
do that consider Eq.4.7, the expression for stress tensor. For hard disk potentials the derivative 
d(f)/dr'^^ becomes a Dirac delta function and the expression for stress can be recast into[102] 



The presence of the Dirac delta function 5{r'^^ — d) in the above expression requires that the terms 
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Figure 4.4: P^i as a function of MC steps. P^i has been multiplied by 10^ and MC steps has been 
multiplied by 10~^ for clarity. The data have been collected for r) = .7029 and pVo = 1. Panel (a) is 
for system size of 2150 particles whereas (b) for 21488 particles. Within 10^ MC steps all fluctuations 
die out. Clearly, the dependence of the dislocation number density on the system sizes and the Monte 
Carlo errors are negligible. To calculate dislocation fugacity we use averaging of data between 5 x 10^ 
to 10^ MC steps. 
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Figure 4.5: Plot of PcPa^x vs. S/d at a strain value e^x = -02 for packing fraction 77 = .7029 and 
potential strength Vq = 1. A second order polynomial fit (solid line) utilizing the error bars to assign 
weights to each data points gives lim^^o Pd'^o'xx = —6.21 with an error within .08%. 

under the summation contribute, strictly, when two hard disks touch each other i.e. r*^ = r = a. 
In practice, we implement this, by adding the terms under summation when each pair of hard 
disks come within a small separation r = a + 5. We then find fiax^d^ as function of 5 and fit 
the curve to a second order polynomial. Extrapolating to the 5^0 limit obtains the value of a 
given component of stress tensor at each strain value eA,/[102](Figs 4.5,4.7). 
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Figure 4.6: A typical stress-strain curve used to obtain the Young's modulus from a linear fit (solid 
line). The graph is plotted at ?? = .7029, Vq = 1.0. The fitted Young's modulus I3K(P = 54.5 with 
an error within 2.9%. The error bars in stress are less than .2% and much smaller than the point sizes 
plotted in this graph. 




Figure 4.7: Plot of (3(faxy vs. 5/d at strain value exy = .079 at the packing fraction r/ = .7029 
and potential strength Vq = 1. A second order polynomial fit (solid line) gives \ims-^o f3(Paxy = 1.033 
with an error within .5%. 



For completeness, now we show how we calculate the two relevant elastic moduli from the 
stresses : a^x at a given longitudinal strain exx (Fig. 4.5) and a^y for a shear strain e^y (Fig. 4.7). 
All the data points are from our MC simulations averaged between 10,000-20,000 MC steps. 
Increasing the number of configurations does not change the values significantly. The total errors 
arising from the MC simulations and the fit for a typical calculation of stresses is less than a 
percent. We thus calculate the stress at each value of strain and from the slopes of stress-strain 
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Figure 4.8: A typical stress-strain curve used to obtain shear modulus from a linear fit (solid line). 
The graph is plotted at 77 = .7029, Vq = 1.0. The fitted shear modulus /?//d^ = 13.5 with an error 
within .9%. The error bars in stress are less than .5% and much smaller than the point sizes plotted in 
this graph. 

curves find out the bare Young's modulus PKd'^ (Fig. 4.6) and shear modulus /3//d^(Fig. 4.8). 
We impose an elongational strain in x- direction which is parallel to the direction of potential 
minima to obtain PKd^. Imposition of a shear strain in the same direction gives us /3//d^. Any 
strain that forces the system to ride potential hills will give rise to massive displacement modes 
which do not contribute to elastic theory. Our results for the stress strain curves for obtaining 
(3Kd'^ and P/id'^ are shown in Figs 4.6,4.8 respectively. Note that the errors for the calculation 
of the elastic constants arise solely from the fitting of the stress- strain curves. These can be 
made as small as possible by increasing the number of strain values at which the stresses are 
calculated. The values of the stress are also free from any residual finite size effects which we 
checked by simulating systems of sizes 10 x 10 to 136 x 158. From these elastic moduli we get 
the 'bare' K,_i.y (and hence Xq = 7rKxy — 2, see section 4.1). Our final estimate for the 'bare' K^y 
are also correct to within a percent. 

In Ref.[95] it is argued that the elastic constant PKd'^ remains more or less independent of 
amplitude of the laser potential (3Vo, while the shear modulus decreases linearly with increasing 
PVo for large pVo- in figures 4.9 and 4.10 we have plotted the values of PKd"^ and Pfid"^ 
respectively as a function of pVo- it is apparent that the expectations of Ref.[95] are borne out 
by our data, incidentally, the behaviour of PKd'^ and Pfid"^ with increasing (3Vo offers an intuitive 
interpretation of the RLiF transition which we offer below. 

in Fig. 4.11 we have plotted x'q and i/q the bare values of x' and y' for various potential 
strengths pVo at packing fraction rj = .7029 along with the separatrices for the linearized and 
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Figure 4.9: Young's modulus PKd^ as a function of inverse laser potential (/3Vb)~^- Various 
symbols denote different densities - o denotes rj = .7029, A denotes 77 = .705 and □ denotes 77 = .7. 
The data for Figs 4.9 and 4.10 were obtained from a separate run with a slightly higher error than that 
in Figs. 4.6 and 4.8. 
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Figure 4.10: Shear modulus /3/id^ as a function of inverse laser potential {f3Vo)~^. Various symbols 
denote different densities - o denotes 77 = .7029, A denotes rj = .705 and □ denotes r] = .7. The 
dotted line is a linear fit of the form /?)ud^ = a/pVo + 6 in the large /9Vb limit[95]. 



the non-linear flow equations (Eq. 4.2). The line of initial conditions is seen to cross the non- 
linear separatrix twice (signifying re-entrant behaviour) while crossing the corresponding linearized 
separatrix only once at high potential strengths. For small pVo the freezing transition is seen 
to be driven mainly by the decrease of y' (the dislocation density) since PKd^ and Pfid"^ are 
virtually constant. For large /?Vo, the shear modulus Pfid"^ vanishes and this results in the second 
point of intersection with the separatrix (remelting). The phase diagram (Fig. 4.13) is obtained 
by computing the points at which the line of initial conditions cut the non-linear separatrix using 
a simple interpolation scheme. It is interesting to note that within a linear theory the KT flow 
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Figure 4.11: The initial values of and y' obtained from the elastic moduli and dislocation 

probability at r/ = .7029 plotted in x' — y' plane. The line connecting the points is a guide to eye. The 
arrow shows the direction of increase in I3Vq{= .01, .04, .1, .4, 1, 4, 10, 40, 100). The dotted line denotes 
the separatrix [y' = x') incorporating only the leading order term in KT flow equations. The solid curve 
is the separatrix when next to leading order terms are included. In Z — oo limit any initial value below 
the separatrix flows to y' = line whereas that above the separatrix flows to y' oo. The intersection 
points of the line of initial values with the separatrix gives the phase transition points. The plot shows 
a freezing transition at /3Vb = -1 followed by a melting at f3Vo = 30. 

equations fail to predict a RLIF transition. Performing the same calculation for different packing 
fractions r] we find out the whole phase diagram of RLIF in the r]- (3Vo plane. 

Small, residual numerical errors in x'q and i/q translate into errors in the location of the phase 
transition points. These are calculated as follows. The quantity PK.j.y varies linearly with r] at all 
potential strengths. Therefore the numerical error in r] is proportional to the error in PK^y (see 
Fig. 4.12). The error in is neglected ^. The final error estimates are shown (as vertical error 
bars) in our results for the phase diagram of hard disks in an external potential in Fig. 4. 13. 

Comparing with previous computations[88, 89] of the phase diagram for this system (also shown 
in Fig. 4.13) we find that, within error- bars, our results agree at all values of r] and PV^ with 
the results of Strepp et. al. [89]. In numerical details, they, however, disagree with the results 
of C. Das et. a/.[88], though even these results show RLIF and are in qualitative agreement with 
ours. This validates both our method and the quantitative predictions of Ref. [94, 95]. 

^Thc error in t/q has contributions from both PdiS and from K^y while the former is practically zero, the 
contribution from the latter is neglected in this work. 
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Figure 4.12: For hard disk system K^y = PK^y varies linearly with rj. Data plotted at Vq = 1. 

The solid line is a linear fit to the form f{x) = a + bx with a = —7.37 and b = 11.76. At each Vq the 
error in K^y determines the error in 77: Srj/r] = |1 + a / r]b\{6Kxy / K^y) ■ 

The effect of higher order terms in determining non-universal quantities has been pointed out 
earlier[68] for the planar rotor model but in the present case their inclusion appears to be crucial. 
Nevertheless, we expect our procedure to break down in the I3Vq limit where effects due 
to the cross-over from a KT to a KTHNY[62, 63] transition at I3Vq = become significant. 
Indeed, as is evident from Fig. 5.12 for pVo < 1 the dislocation probabilities of both type I and 
type II dislocations are similar ^ and the assumptions of FNR theory and our process (which 
involves only type I dislocations) need not be valid at small potential strengths. This fact is also 
supported by Ref.[89] where it was shown that though at pVo = 1000 the scaling of susceptibility 
and order parameter cumulants gave good data collapse with values of critical exponents close 
to FNR predictions, at /3Vo = .5, on the other hand, ordinary critical scaling gave better data 
collapse than the KT scaling form, perhaps due to the above mentioned crossover effects. In the 
asymptotic limit of PVq — > 00 the system freezes above 77 = .706 which was determined from a 
separate simulation in that limit. This number is very close to the earlier value 77 ~ .71 quoted in 
Ref.[89]. As expected, the freezing density in the pVo — > 00 limit is lower than the value without 
the periodic potential i.e. rj ~ .715. 

Before we go on to discuss other systems, we discuss the reasons behind the particular choice 

^In analysing Fig. 3 we must keep in mind that we can calculate from our simulations only the probability 
of formation of a disclination quartet. While we can, perhaps, safely assume that if type I dislocations 

arc involved, they will scpcratc out with unit probability, the same can not be said of type II dislocations. 
This means that the probability of type II dislocations could be much lower than what Fig. 3 suggests. 
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1000 

Figure 4.13: The phase diagram of the hard disk system in the presence of a 1-d, commensurate, 
periodic potential in the packing fraction (ry) - potential strength (/3Fo) plane. The points denoted by □ 
correspond to our RG calculation using the techniques described in this chapter. The points denoted by 
0[89] and *[88] are taken from earlier simulations. The vertical bars denote estimate of error. Our data 
clearly matches with Ref[7].The horizontal line at = .706 denotes the calculated asymptotic phase 
transition point at pVo = oo. 

of r-min that we made throughout this work. After a disclination quartet is formed, they get 
separated out and the easy direction of separation is the glide direction which is parallel to the 
Burger's vector. In Fig. 4. 14 we show four steps of separation of such a dislocation pair of type 
I. It is clear that, it is possible to give individual identification to a dislocation only when the 
Burger's vector separation within a pair is > 2ao (Fig. 4. 14) i.e. r-min = 2ao. For r > 2ao 
Burger's loops can be drawn around each 5 — 7 disclination pair (Fig. 4. 14) giving rise to a non- 
zero Burger's vector. After motivating r^m = 2ao we show, in Fig. 4. 15, the three sets of initial 
values corresponding to r^jn = ao, 2ao, 3ao along with the non-linear separatrix at 77 = .7029 
of hard disk system. It is clear from the figure that r^m = o-o predicts the system to be in the 
solid phase for any arbitrarily small amount of external potential and to melt at larger /dVo. This 
behaviour contradicts physical expectation that the melting density at pVo — has to be larger 
than that at pVo — 00. On the other hand, while rmin — 3ao does not produce any unphysical 
prediction, it shrinks the region of re-entrance in the pVo direction. It is therefore satisfying to 
note that r^j„ = 2ao, the minimum possible value for the separation between members of a 
dislocation- antidislocation pair which allows unambiguous identification also produces physically 
meaningful results for the phase diagram in closest agreement with earlier simulation data. 

It is possible to find out phase diagrams of any 2-d system in presence of external modulating 
potential commensurate with the density of the system in a similar fashion. We illustrate this by 
calculating similar phase diagrams for two other systems, viz. soft disks and the DLVO system. 
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Figure 4.14: The figures a - d which we have drawn using the applet " voroglide" [111] show four 
steps of separation of a type I dislocation pair, from a separation of uq to Aqq. The shaded regions 
show the 5 — 7 disclination pairs constituting the dislocations. Burger's circuits are shown in a - c. 
Note that for separations > 2ao separate Burger's circuits around each disclination pair give rise to non- 
zero Burger's vectors, giving the dislocations their individual identity. This shows that the minimum 
meaningful separation between dislocation cores Vmin = 2ao. 

Soft disks: Soft disks interact via the potential : 

'^(^) = :^ 

where r denotes the separation between particles. In simulations, the cutoff distance is chosen to 
be Tc = 2 above which the particles are assumed to be non- interacting. Apart from the external 
potential strength /3Vo the relevant thermodynamic quantity is the number density p — N/L^Ly. 
To obtain 'bare' elastic moduli from restricted simulations the stress is calculated using Eq.4.7. As 
this expression does not involve any Dirac delta functions (unlike hard disks), we do not require 
any fitting and extrapolation to obtain the stresses and the errors are purely due to random 
statistical fluctuations in our MC simulations. The elastic moduli are again found from stress- 
strain curves like Figs.4.6,4.8. The dislocation fugacity of type I is calculated from rejection ratio 
of dislocation generating moves. All these, at a given p value generate the initial conditions x'q 
and i/o in RG flow diagrams. Again, the crossing of these initial conditions with the separatrix 
found from Eq.4.2 gives the phase transition points. The phase diagram is plotted and compared 
with phase diagram from earlier simulations[91, 92] in Fig. 4. 16. The error bar in p is found from 
the error in K^y, as K^y varies linearly with p, through the relation 6p/p = \\ + a/ ph\{5Kxy/ K^y) 
where a and h are obtained from a linear (a + hx) fit of the K^y vs. p curve, at any given BVq. 
The phase diagram (Fig. 4.16) again clearly shows re-entrance (RLIF). This is in qualitative 
agreement with earlier simulations[91, 92] (see Fig. 4. 16). 

DLVO: For charge stabilized colloids the inter-particle potential that operates is approximately 
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Figure 4.15: Similar to Fig. 4. 11. The initial conditions x'q and are plotted as a function of PVq. 
The different data sets are created for different values of rmin- The synnbols mean the following : 
denotes data for r^m = oq, □ denotes that for rmin = 2ao and A denote data for rmin = 3ao. The 
dotted line denotes the non- linear separatrix. 
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Figure 4.16: Phase diagram for soft disks: □ denote our calculation, indicate earlier simulation 
data [91, 92]. The vertical lines are the error- bars. 



given by the DLVO potential [103, 104]: 

, . {Z*ey f exp{.hK(i)\^ exp{—Kr) 
~ A-Ke^er V 1 + .Sred ) r 

where r is the separation between two particles, d is the diameter of the colloids, k is the inverse 
Debye screening length, Z* is the amount of effective surface charge and is the dielectric 
constant of the water in which the colloids are floating. In order to remain close to experimental 
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Figure 4.17: Phase diagram for particles interacting via the DLVO potential. □ denote our calcu- 
lation, show the earlier simulation data[90]. The vertical lines are the error- bars. Error bars in our 
calculation being smaller than the symbol size are not shown. 



situations and to be able to compare our phase diagram with the simulations of Strepp et. al.[89] 
we use T = 293. 15i^', d= 1.07/im, Z* = 7800, = 78. In experiments, the dimensionless 
inverse Debye screening length nug can be varied either by changing k through the change 
in counter-ion concentration or by changing ag by varying density[112]. In our restricted MC 
simulations we vary k, keeping the density fixed at 0.18/xm~^ by fixing the lattice parameter of 
the in ital configuration of ideal triangular lattice at = 2.52578/.tm. Further, we use a cut- 
off radius Tc such that, 0(r > Vc) — where Tc is found from the condition /30(rc) — .001. 
We find out phase transition points (in nag) at different external potential strengths PVq in 
the same fashion as described earlier. The phase diagram in KUg — pVo plane is shown in 
Fig. 4.17. To obtain error bars in this case we note that Kxy varies linearly with KUg and 
therefore the error in K^y is proportional to the error in Kttg (Fig. 4. 17) through the relation 
5{K,as)/{K,as) = |1 + a/bK,as\{5Kxy/ Kxy). The quantities a and b are found from fitting K^y to 
a linear form of nag, at any given PVq. 

Though there is a quantitative mismatch between our data and that of Strepp et. a/. [90], our 
data shows a clear region in KUg (between 15.1 and 15.2) where we obtain re-entrance (RLIF). 
This is in contrast to the simulated phase diagram of C. Das et a/.[88], where they observe 
absence of re-entrance at high field strengths. We do not plot their data as the parameters these 
authors used are not the same as the ones used in Fig. 4. 17. 

It is interesting to note that, with increase in the range of the two- body interaction potentials 
the depths of re-entrance (in r], p or nag) decreases. This is again in agreement with the 
understanding that, the re-entrant melting comes about due to decoupling of the 1-d trapped 
layers of particles that reduces the effective dimensionality thereby increasing fluctuations. With 
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an increase in range of the interacting potentials this decoupling gets more and more suppressed, 
thereby reducing the region of re-entrance. 

One aspect of our study which stands out is the exceptionally better agreement of our results 
with previous simulations for hard disks as opposed to systems with soft potentials like the soft 
disks and the DLVO. This could, in principle, be due either (a) to the failure of the RG equations 
used by us or some other assumptions in our calculations (b) or to unaccounted finite size effects 
in earlier simulations. While it is difficult to estimate the effect of (a) since RG equations to 
higher orders in y are unknown at present, we may be able to motivate an estimation for (b). In 
order to explain the discrepancy in the positions of the phase boundaries, we need to go into some 
details of how the phase diagrams were obtained in the earlier simulations. In these simulations 
[89-92] the phase boundaries were obtained from the crossing of the order parameter cumulants 
[5, 97] for various coarse graining sizes. The system sizes simulated in these studies are the same 
{N — 1024). However, the range of interaction differs. To obtain an objective measure we define 
the range of the potentials ^ as that at which the interaction potential is only 1% of its value 
at the lattice parameter. In units of lattice parameter, we obtain, for soft disks ^ = 1.47 and 
for the DLVO potential ^ = 1.29 at typical screening of — 15. By definition, for hard disks 
^ = 1. The particles within the range of the potential are highly correlated and we calculate 
the number Ncorr of such independent bare uncorrelated particles within the full system size. 
Ncorr takes the values N^orr = 1024, 473.88, 615.35 for hard disks, soft disks and the DLVO 
potential respectively. Since the effective system sizes are smaller for the soft potentials, finite 
size effects are expected to be larger. In this connection, it is of interest to note that in the 
same publications [89-92] a systematic finite size analysis showed that the phase diagrams shift 
towards higher (lower) density (kappa) for hard and soft disks (DLVO). A look at Fig. 4. 16 and 
4.17 should convince the reader that such a shift would actually make the agreement with our 
results better. We emphasize here that our present restricted simulations are virtually free of 
finite size effects since the system does not undergo any phase transition. 

4.3 Conclusion 

We have presented a complete numerical renormalization group scheme to calculate phase dia- 
grams for 2-d systems under a commensurate modulating potential. We have used FNR theory 
along with this scheme to calculate phase diagrams for three different systems, namely, the hard 
disks, the DLVO and the soft disks. In all the cases we have found laser induced freezing followed 
by a re-entrant laser induced melting. We show that the re-entrance behavior is built into the 
'bare' quantities themselves. We find extremely good agreement with earlier simulation results. 
In particular the phase diagram for hard disk comes out to be exactly the same as found from 
one set of earlier simulations[89].To obtain the correct phase diagram, however, flow equations 
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need to be correct at least upto next to leading order terms in the dislocation fugacity. Our 
results, especially for small potential strengths, is particularly sensitive to these terms. Cross-over 
effects from zero potential KTHNY melting transition are also substantial at small values of the 
potential. 

In this chapter we have studied the phenomena of RLIF, that comes about due to a confining 
potential which is constant in one direction and modulating in the other. In next chapter we shall 
study the effect of another kind of confinement. We shall confine a two dimensional solid in a 
narrow but long channel and will find out its properties, phases, strain induced failure and phase 
diagram. 
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5 Confined Solid: Phases and Failure 



One afternoon the boys grew enthusiastic over the flying carpet that went swift:ly 
by the laboratory at the window level ... - G. G. Mdrquez 



In the last chapter we studied phase transitions in a two dimensional solid driven by a periodic 
confining potential which caused a dinnensional crossover from two to one dimension as the 
amplitude of the periodic potential was increased. In this chapter we study the effect of a different 
kind of external potential which forces a system of "hard-disk" atoms to remain in-between one 
and two dimensions. Specifically, we consider here the mechanical and thermodynamic properties 
of a narrow strip of crystalline solid with one dimension much longer than the other other i.e. 
a quasi one dimensional (QID) system. In the short dimension, the solid is confined by hard, 
featureless walls. 

We shall show in this chapter, that such a QID solid strip has rather anomalous properties, which 
are quite different from bulk one, two or three dimensional systems. The QID solid is shown to 
have a non zero Young's modulus which offers resistance to tensile deformations and approximate 
two dimensional hexagonal crystalline order. On the other hand, the shear modulus of the system 
is vanishingly small. Large wavelength displacement fluctuations are seen to destabilize crystalline 
order beyond a certain length scale at low densities. At high densities these fluctuations appear 
to be kinetically suppressed. The failure properties of this quasi solid is also rather interesting. 
In the constant extension (Helmholtz) ensemble, the initial rise in the tensile stress with tensile 
loading is interrupted at a limiting value of strain and on further extension the stress rapidly falls 
to zero accompanied by a reduction in the number of solid layers parallel to the hard wall by one. 
However, this failure is reversible and the system completely recovers the initial structure once 
the strain is reduced. The critical strain for failure by this novel mechanism, for small channel 
widths, decreases with increase in channel width so that thinner strips are more resistant to 
failure. We have used an idealized model solid to illustrate these phenomena. Our model solid 
has particles (disks) which interact among themselves only through excluded volume or "hard" 
repulsion. We have reasons to believe, however, that for the questions dealt with in this chapter, 
the detailed nature of the inter particle interactions are relatively irrelevant and system behaviour 
is largely determined by the nature of confinement and the constraints. Our results may be 
directly verified in experiments on sterically stabilized "hard sphere" colloids[96] confined in glass 
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channels and may also be relevant for similarly confined atomic systems interacting with more 
complex potentials. We have also speculated on applications of this reversible failure as accurate 
strain transducers or strain induced electrical or thermal switching devices. In the next chapter 
we shall study thermal transport across this QID solid, especially with respect to the effects of 
reversible failure on the transport coefficients. 

Studies of small assemblages of molecules with one or more dimensions comparable to a few 
atomic spacings are significant in the context of nano-technology[113, 114]. Designing nano- 
sized machines requires a knowledge of the mechanical behavior of systems up to atomic scales, 
where, a priori, there is no reason for our ideas, derived from macroscopic continuum elasticity 
theory, to be valid[115]. Small systems often show entirely new behaviour if hard constraints are 
imposed leading to confinement in one or more directions. Consider, for example, the rich phase 
behaviour of quasi two-dimensional colloidal solids [116-121] confined between two glass plates 
showing square, triangular and "buckled" crystalline phases and a, recently observed, re-entrant 
surface melting transition[122] of colloidal hard spheres not observed in the bulk[71, 99, 105-107]. 

Recent studies on various confined QID systems have shown many different structures depending 
on the range of interactions and commensurability of the natural length scale of the system with 
the length scale of confinement [123, 124]. These structures play crucial role in determining the 
local dynamical properties like asymmetric diffusion, viscosity etc. and phase behaviour [125]. 
A study by G. Piacente et. al. [123] on confined charged particles interacting via the screened 
Coulomb potential and confined in one direction by a parabolic potential showed many zero 
temperature layering transitions, i.e. change in the number of layers with a change in density or 
range of interaction. This also showed regions where these transitions were reentrant [123]. Apart 
from the transition from one to two layers, all these layering transitions were shown to be first 
order [123]. At high temperature this classical Wigner crystal melts and the melting temperature 
shows oscillations as a function of density [123]. Such oscillations are charecteristic of confined 
systems, arising out of commensurability. Confined crystals always align one of the lattice planes 
along the direction of confinement [123, 124] and confining walls generate elongational asymmetry 
in the local density profile along the walls even for the slightest incommensuration. A study by R. 
Haghgooie et. al. on a system of purely repulsive dipoles confined in QID hard channel showed 
layering transitions mediated via an order- disorder transition near the centre of the channel [124]. 
All the structural properties showed oscillations as the channel width increased [124]. Wall 
induced layering was also observed in a dusty-plasma study by L.-W. Teng et. al. [126] and in 
shell structures in circular confinements [127-131]. A similar layering transition, in which the 
number of smectic layers in a confined liquid changes discretely as the wall-to-wall separation is 
increased, was noted by P. G. deGennes [132] and J. Gao et. al. [133]. J. Gao et. al.'s study also 
revealed the relations between the local layering structures and dynamical quantity like diffusion 
constant [133]. deGennes and Gao et. al.'s work also calculated the wall to wall force due to 
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this smectic layering and its commensurability with channel width. For long ranged interactions 
extreme localization of wall particles has been observed in studies of V. M. Bedanov et. al. [127] 
and R. Haghgooie et. al. [124]. In an earlier experiment on confined steel balls in QID vibrated 
to simulate the effect of temperature, layering transition, phase coexistance and melting was 
observed by P. Pieranski et. al. [134]. Confinement, in general, can lead to new behaviours 
in various systems and processes. Schmidt and Lowen [120] studied the phase behaviour of a 
collection of hard spheres confined within a two dimensional slit defined by two parallel hard 
plates. Plate separations such that upto two layers of solid are accomodated were considered by 
these authors. Recently Fortini and Dijkstra studied the same system for separations of one to 
five hard sphere diameters to find a rich phase diagram consisting of a dazzling array of upto 
26 distinct crystal structures[121]. Similar studies had been carried out long back in 1983 by 
Pieranski and his group [116] to find out many of the structures identified by Fortini et. al. [121]. 
In biology specific structures of proteins are required for their specific functioning which keeps 
a cell living. Double barreled chaperonins capture protiens into them to fold them to specific 
structures. Recently, an array of genetically engineered chaperonin templates have been used by 
McMillan et. al. [135] to produce ordered nano- particle arrays. 

This chapter is organized as follows. In the next section, we shall introduce the model confined 
solid and discuss the geometry and basic definitions of various structural and thermodynamic 
parameters. We shall then introduce the various possible structures and phases with their basic 
characteristics. This will be followed by the results of computer simulations, in the constant 
NAT (number, area, temperature) ensemble, exploring the deformation and failure properties 
of this system and the relation of the various structures described in the previous section to 
one another. In the fifth section we shall try to understand the various transitions seen in the 
computer simulations within simple mean field free volume and density functional approaches. In 
section VI we shall show the effect of fluctuations and its role in destruction of long range order 
in this system. We conclude the chapter in section VII. 

5.1 The Model System 

The bulk system of hard disks where particles i and j, in two dimensions, interact with the 
potential Vij = for \rij\ > d and Vij = oo for \rij\ < d, where d is the hard disk diameter and 
Yij = Tj — Ti the relative position vector of the particles, has been extensively[71, 99, 105-107] 
studied. Apart from being easily accessible to theoretical treatment[109], experimental systems 
with nearly "hard" interactions viz. sterically stabilized colloids[96] are available. The hard 
disk free energy is entirely entropic in origin and the only thermodynamically relevant variable 
is the number density p = N/V or the packing fraction r] = (7r/A)pd'^. Accurate computer 
simulations[99] of hard disks show that for rj > rjf = .719 the system exists as a triangular 



68 



Confined Solid: Phases and Failure 




Figure 5.1: The confined solid is shown along with the centered rectangular (CR) unit cell. For 

an unstrained triangular lattice Ux = ao and ay = \/3ao. Gi, G2 and G3 denote the directions of 
the three reciprocal lattice vectors (RLV). The third reciprocal lattice direction G3 is equivalent to the 
direction G2, even in presence of the walls. 



lattice which melts below 77^ = .706. Elastic constants of bulk hard disks have been calculated 
in simulations[71, 136]. The surface free energy of the hard disk system in contact with a hard 
wall has also been obtained[137] taking care that the dimensions of the system are compatible 
with a strain-free triangular lattice. 

Consider a narrow channel in two dimensions of width Ly defined by hard walls at y = and 
Ly (Kraii(y) — for d/2 < y < Ly — d/2 and = 00 otherwise) and length with ^ Ly. 
Periodic boundary conditions are assumed in the x direction(Fig.5.1). In order that the channel 
may accommodate ni layers of a homogeneous, triangular lattice with lattice parameter oq of 
hard disks of diameter d, (Fig. 5. 7) one needs, 



-{ni - l)ao + d 



(5.1) 



For a system of constant number of particles and Ly, oq is decided by the packing fraction rj 
alone. Note that L^ = n^^ao = Nao/rii, and oq is given by p = N/L^Ly. This gives 



=^ + ,/4 + 2v^(l- ^)i 



V3(l - 



(5.2) 



Defining xiv^^y) = 1 + "^{Ly — d)/^/3ao, the above condition reads x = integer = rii and 
violation of Eqn.(5.1) implies a rectangular strain away from the reference triangular lattice of 
ni layers. The lattice parameters of a centered rectangular (CR) unit cell are and ay (Fig. 
5.1). In general, for a CR lattice with given Ly we have, ay = 2{Ly — d)/{ni — 1) and, ignoring 
vacancies, a^ = 2/ pay. 
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Figure 5.2: A plot of internal strain as a function of external strain e. The jumps in 
corresponds to half-integral values of x- 



Calculation of the deformation strain needs some care at this stage. Using the initial triangular 
solid (packing fraction -qo) as reference, the "external" strain associated with changing L^, while 
keeping and Ly fixed, is e = (L^ — = {jiq — r])/r] where t] is the packing fraction of 

the deformed solid and r] = Nud^ / ALyL.^.. Internally, the solid is, however, free to adjust rii to 
decrease its energy (strain). Therefore, one needs to calculate strains with respect to a reference, 
distortion-free, triangular lattice at r]. Using the definition Sd = Sxx —^yy = i^x — do) /do — (dy — 
V3ao)/\/3ao = a^/ao — ay/^/Sao and the expressions ax = '2/ pay, ciy = '^{Ly — d)/{ni — 1), 
ao = 2(Lj, — d)/\/3{x — 1) we obtain, 



ni-l _ X - 1 
X - 1 ni-l' 



(5.3) 



where the number of layers ni is the nearest integer to x so that has a discontinuity at half - 
integral values of x- For large Ly this discontinuity and Sa itself vanishes as '\-/ Ly for all rj. This 
"internal" strain is related non-linearly to e and may remain small even if e is large (Fig. 5. 2). 
Note that any pair of variables 77 and Ly (or alternately e and x) uniquely fixes the state of the 
system. 



5.2 Structures and Phases 



The different possibilities of structures and phases along with their typical structure factors are 
presented in this section. 
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Figure 5.3: Solid: Left panel shows an picture of 10^ overlapped configurations of a high densitiy 
[t] = .85) solid phase with wall to wall separation commensurate with the density. The color code is such 
that red means high local density and blue means low density. The right panel shows the corresponding 
structure factor which shows the typical pattern for a two dimensional hexagonal solid. 



If the separation between the hard walls is kept commensurate such that x = ^i, some integer 
number of layers then the equilibrium phase is a perfect two dimensional triangular solid (Fig. 5. 3). 
The solid shows a diffraction pattern which is typical of a two dimensional trangular crystal. We 
show later that appearances can be deceptive, however. This triangular "solid" is shown to have 
zero shear modulus which would mean that it can flow without resistance along the length of the 
channel like a liquid. Stretching the solid strip lengthwise, on the other hand, costs energy and is 
resisted. The strength of the diffraction peaks decreases rapidly with the order of the diffraction. 
In strictly two dimensions this is governed by a nonuniversal exponent dependent on the elastic 
constants [62]. In QlD this decay should be faster but larger system sizes and averaging over a 
large number of configurations would be required to observe this decay since contraints placed by 
the hard walls makes the system slow to equilibriate at high densities. We return to this question 
in section VI. 

As discussed in the earlier section, even a small incommensuration due to the confining walls in 
QlD immediately introduces elongational asymmetry to local density profiles along the confining 
directions. As a result of this, a nonzero elongational stress is induced in the system. This 
causes the two diffraction spots corresponding to planes parallel to the hard walls to strengthen 
at the cost of the other four spots in the smallest reciprocal lattice set. This increases the one 
dimensional character of the system even further. 

A little extra space introduced between the walls starting from a high density solid phase gives 
rise to buckling instability in y- direction and the system breaks into triangular solid regions along 
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Figure 5.4: Buckled phase: A small incommensuration of wall to wall separation introduced via 
a small increase in it from a separation that is commensurate with a high density triangular solid at 
T] = 0.89 gives rise to this phase. Increase in channel width reduces the density to rj = 0.85. The left 
panel shows the picture corresponding to lO'^ overlapped equilibrium configurations. The color code for 
the local denisties are same as before. Note different portions of triangular solid separated in x- direction 
are displaced along y- direction to span the extra space introduced between the walls. Lines are drawn 
to identify this shift in triangular regions. The right panel shows the corresponding structure factor 
which has the peak in Gi direction diminished. Some extra weak peaks corresponding to superlattice 
reflections appear at lower values of the wavenumber. 



the X- direction (Fig. 5. 4). Each of these regions fluctuate with respect to the other in y- direction 
giving the impression of a buckling wave travelling along the length of the solid. In conformity 
with the two dimensional analog [118-120] we call this the buckled solid and it interpolates 
continuously from x = ri; to ri; ± 1 layers. This phase can also occur due to the introduction of 
a compressional strain in x-direction keeping Ly fixed. We do not observe the buckled solid at 
low densities close to the freezing transition. Extreme incommensuration at such densities lead 
to creation of bands of the smectic phase within a solid eventually causing the solid to melt (see 
next section). The diffraction pattern shows a considerable weakening of the spots corresponding 
to planes parallel to the walls, together with an extra spot at smaller wavenumber corresponding 
to the buckled superlattice. The diffraction pattern is therefore almost complementary to that of 
the smectic phase to be discussed below. 

At low enough densities or high enough incommensuration the elongated density profiles in the 
lattice planes parallel to the walls can overlap to give rise to a smectic phase (Fig. 5. 5) in which 
local denity peaks are completely smeared out in x- direction but are clearly seperated in y- 
direction giving rise to a solid like order in that direction and making the system liquid like in x- 
direction. The diffraction pattern shows two spots which is typical corresponding to the symmetry 
of a smectic phase. 
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Figure 5.5: Smectic: The left panel shows the picture of 10^ equilibrium configurations. The color 
code for local density is same as before. The middle panel shows the density modulation in y- direction 
in a typical smectic phase which is solid like in a direction perpendicular to the walls (Gi direction) 
and liquid like in the other direction. A fact coroborated by solid- like peak in Gi direction in structure 
factor plotted in right hand panel. This is at a packing fraction t] = 0.73 obtained the straining a 
triangular lattice at rj = 0.85 in x- direction. 

At further lower densities the relative Lindeman parameter, a quantity which measures the rel- 
ative displacement fluctuations between neighbours (will be defined in following section) diverges 
and the structure factor shows a ring- like feature typical of a liquid appears together with the 
smectic like peaks in Gi direction. This is a modulated liquid (Fig. 5. 6). The density modulation 
decays away from the walls and in channels with larger widths, the density profile in the middle 
of the channel becomes uniform. 

5.3 Mechanical Properties and Failure 

A bulk solid, strained beyond it's critical limit, fails by the nucleation and growth of cracks[138- 
141]. The interaction of dislocations or zones of plastic deformation [141, 142] with the growing 
crack tip determines the failure mechanism viz. either ductile or brittle fracture. Studies of the 
fracture of single-walled carbon nanotubes[143, 144](SWCNT) also show failure driven by bond- 
breaking which produces nano cracks which run along the tube circumference leading to brittle 
fracture. Thin nano-wires of Ni are known[145, 146], on the other hand, to show ductile failure 
with extensive plastic flow and amorphization. 

In this section we shall present our results for the mechanical behaviour and failure of the QID 
solid under tension. As mentioned earlier, we show that the QlD solid behaves anomalously, 
showing a reversible plastic deformation and failure in the constant extension ensemble. The 
failure occurs by the nucleation and growth of smectic regions which occur as distinct bands 
spanning the width of the solid strip. 



We study the effects of strain on the hard disk triangular solid at fixed Ly large enough to 
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FiGURE 5.6: Modulated liquid: The upper panel shows the picture of 10^ equilibrium configurations. 
Local denisties are coded in the same color code as before. The lower left hand panel shows the density 
modulation in y- direction which is like the smectic phase but the modulation dies out at the centre. 
The structure factor in right hand panel shows a ring structure which is a typical signature of liquid 
superimposed on smectic- like strong peaks in Gi direction. This is the signature of a modulated liquid. 
This behaviour is found at a packing fraction rj = 0.6 obtained the straining a triangular lattice at 
rj = 0.85 in x- direction. 



accommodate a small number of layers ~ 9 — 25. We monitor the Lindemann parameter 

/ =< {u\ - u'jf > /al+ < {uy^ - uy,f > /al 

where the angular brackets denote averages over configurations, i and j are nearest neighbors 
and u°'i is the a-th component of the displacement of particle i from it's mean position. The 
parameter / diverges at the melting transition [147]. We also measure the structure factor 



Pg 



1 ^ 
— J2 exp(-zG.r,fe) 



j,k=i 



for G = ±Gi(?7), the reciprocal lattice vector (RLV) corresponding to the set of close-packed 
lattice planes of the CR lattice perpendicular to the wall, and ±G2{f]) the four equivalent RLVs 
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Figure 5.7: Results of NVT ensemble Monte Carlo (MC) simulations of iV = x n^, = 65 x 10 
hard disks confined between two parallel hard walls separated by a distance Ly = 9d. For each 77, 
the system was equilibrated over 10^ MC steps (MCS) and data averaged over a further 10^ MCS. At 
rj = 0.85 we have a strain free triangular lattice. Plots show the structure factors pGi,i = l(+))2(o) 
for RLVs Gi(?7), averaged over symmetry related directions, as a function of rj. Also plotted in the same 
graph is the Lindemann parameter !(□). The lines in the figure are a guide to the eye. Inset shows the 
geometry used, the reciprocal lattice vectors (RLVs) Gi and G2 and the CR unit cell. 

for close-packed planes at an angle (= 7r/3 and 2n/3 in the triangular lattice) to the wall (see 
Fig. 5.1). Notice that ±G2(r;) and ±02(77) as shown Fig.5.1 in are equivalent directions. 

Throughout, < PGi 7^ 0, a consequence of the hard wall constraint[137] which manifests 
as an oblate anisotropy of the local density peaks in the solid off from commensuration. As 77 is 
decreased both pci and show a jump at 77 = 77^ where x — X* ~ 7^^ — 1/2 (Fig. 5.8 (inset)). 
For 77 < 77ci we get po^ — with pci signifying a transition from crystalline to smectic 
like order. The Lindemann parameter I remains zero and diverges only below 77 = 7703 (~ 77,^) 
indicating a finite-size- broadened melting of the smectic to a modulated liquid phase. The stress, 
c — (^xx — Cyy versus strain, e, curve is shown in Fig. 5.8. For 77 = ?7o (e = 0) the stress is purely 
hydrostatic with a^x = o'yy as expected. At this point the system is perfectly commensurate with 
channel width and the local density profiles are circularly symmetric. Initially, the stress increases 
linearly, flattening out at the onset of plastic behavior at 77 ~ 77C1. At 77^, with the nucleation of 
smectic bands, a decreases and eventually becomes negative. At r/cj the smectic phase spans 
the entire system and a is minimum. On further decrease in rj towards 1]^.,^ , a approaches from 
below (Fig. 5.8) thus forming a Van der Waals loop. If the strain is reversed by increasing i] back 
to r]o the entire stress-strain curve is traced back with no remnant stress at i] = r]o showing that 
the plastic region is reversible. For the system shown in Figs. 5. 7 and 5.8, we obtained 77^ ~ .77, 
rjc2 ~ -74 and ?7c3 ~ .7. As Ly is increased, 77^ merges with ric.j for ni ~ 25. If instead, and 
Ly are both rescaled to keep x = f^i fixed or periodic boundary conditions are imposed in both x 
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Figure 5.8: A plot of the conjugate stress a versus external strain e obtained from our MC sim- 
ulations of 65 X 10 hard disks initially at r/ = 0.85. Data is obtained by equilibrating at each strain 
value for 2 x 10'* MCS and averaging over a further 3 x 10^ MCS. The stress for the hard disk system 
has been calculated by the standard method[102]. The entire cycle of increasing e(o) and decreasing to 
zero (+) using typical parameters appropriate for an atomic system, corresponds to a real frequency of 
uj lOOKHz. Results do not essentially change for lo = lOKHz — IMHz. Inset shows the variation 
of the critical x* with ni, points: simulation data; line: x* = "^i ~ 1/2- 

and y directions, the transitions in the various quantities occur approximately simultaneously as 
expected in the bulk system. Varying in the range 10 — 1000 produces no essential change in 
results. 

For r]c2 < r] < Vci we observe that the smectic order appears within narrow bands (Fig. 5.9). 
Inside these bands the number of layers is less by one and the system in this range of is in a 
mixed phase. A plot (Fig. 5. 9 (a) and (b)) of xix, t), where we treat x as a space and time (MCS) 
dependent "order parameter" (configuration averaged number of layers over a window in x and 
t), shows bands in which x is less by one compared to the crystalline regions. Once nucleated 
narrow bands coalesce to form wider bands, the dynamics of which is, however, extremely slow. 
The total size of such bands grow as rj is decreased. Calculated diffraction patterns (Fig. 5.9 
(c) and (d)) show that, locally, within a smectic band pci ^ PCa in contrast to the solid region 
where pci ~ Pg-, 7^ 0. 

The total free energy per unit volume of a homogeneous soWd, , which is in contact with a 
hard wall and distorted with a (small) strain is given by. 



(5.4) 
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Figure 5.9: Plot of xix,t) as a function of the x/d at r/ = .76 after time t = (a)5 x 10^ and 

(b) 2 X 10^ MCS for iV = 10^ x 10. Note that ^ = 10 in the solid and = 9 in the smectic regions. 
Arrows show the coalescence of two bands as a function of time. Calculated diffraction patterns for the 

(c) solid and (d) smectic regions, (e) Close up view of a crystal-smectic interface from superimposed 
positions of 10'^ configurations at 77 = .77. The colors code the local density of points from red/dark 
(high) to blue/light (low). Note the misfit dislocation in the inter-facial region. 



where K^{r]) is an elastic constant and J-'^{i]) the free energy of the (undistorted) triangular 
lattice in contact with a hard wall[137] at packing fraction r]. The "fixed neighbor" free volume 
Vf{r], x) may be obtained using straight forward, though rather tedious, geometrical considerations 
so that J^^(r7) = -kBTp\ogVf{ri,0) and K^{ri) = d^J^^{ri,ed)/del\e^=o (see Fig.5.8). Vf is 
expressed in units of d'^. In fixed neighbour free volume theory (FNFVT), we think of a single 
disk moving in a fixed cage formed by taking the average positions of its six nearest neighbor 
disks [ see Fig. (5.10) ]. The free volume available to this central particle is given entirely 
by the lattice separations b (= AQ) and h (See Fig. 5. 10). Note that, b = ao(l + txx) and 
h = -\/3(l + tyy)/2 where ao is lattice parameter of a triangular lattice at any given packing 
fraction 77 and e^^ = ~ 1). ^yy = {x~ ~ !)■ As stated in Sec. 5.1, x is obtained 
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Figure 5.10: In our free-volume theory we assume that the outer six disks are fixed and the central 
disk moves within this cage of fixed particles. The curve in bold line shows the boundary B of the free 
volume. A point on this boundary is denoted by Po{x,y) while the centers of the six fixed disks are 
denoted by Pi{xi,yi) with i = 1,2. ..6. 



from channel width Ly and packing fraction 77. Vf is the area enclosed by the boundary B in 
Fig. 5. 10. In calculating the free volume we have assumed the geometry of the free volume to 
be close to hexagonal (Fig.5.10). For small strains (within 6%) around a triangular lattice this 
assumption is valid. However, for large strains the free volume area goes over to a rhombic shape. 
At these points our theory fails. It is clear that JF^ has minima for all x = n-i. For half integral 
values of x the homogeneous crystal is locally unstable. The FNFVT fails also at these points. 
This same FNFVT free energy for solid is used in Sec.5.5 to calculate the phase diagram of this 
confined system. 

Noting that, x* = — 1/2 (Fig. 5.8 inset), it follows from Eqn. 5.3, the critical strain 
e*^ — {Aril — 5)/{'2ni — 3)(2n; — 2) ~ 1/n; which is supported by our simulation data over the 
range 9 < ni < 14. At these strains the solid generates bands consisting of regions with one less 
atomic layer. Within these bands adjacent local density peaks of the 'atom's overlap in the x 
direction producing a smectic. Indeed, the overlap maybe calculated approximately using simple 
density functional arguments[148] to be A = ^J<~u^/> j = (x — l)/47rA/CoPG7("'/^l) (where 
Co, direct correlation function for a hard disk uniform liquid, is a constant of order unity) which, 
evidently, diverges as (see Sec.5.5). Note that at x = X* the term (x — 1)/ {ni — 1) and 

hence the overlap shows a jump discontinuity even before ~^ 0. In Fig. 5. 11 this overlap term 
has been plotted as a function of strain e imparted on a ten layered solid. Fig. 5. 11 shows that the 
overlap A undergoes a jump increase of around 150% at e ~ 0.1 and reaches its highest value 
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Figure 5.11: Values of pG2. X st^d ni have been calculated from the same simulation which is used 
in plotting Fig.5.7. Co has been set to one in calculating density profile overlap A at each strain value 



at e ~ 0.15 as the smectic phase spans the whole system. At further higher strain this overlap 
decreases with melting of the smectic into fluid phase. This observation further vindicates the 
phase demarcation in Fig. 5. 8. For large Ly, the failure strain e*^ reduces significantly to wash 
out the difference in maxima and minima in free energy, therefore, the minima in J^^ merge to 
produce a smooth free energy surface independent of x and more conventional modes of failure, 
viz. cracks, are expected to become operative. 

For small Ly all regions of the parameter space corresponding to non-integral x s'so globally 
unstable as belied by the loop in the stress-strain curve (Fig. 5. 8). The system should therefore 
break up into regions with ni and n/ — 1 layers for infinitesimal Sd- Such fluctuations are, however, 
kinetically suppressed as we argue below. A superposition of many particle positions near such 
an interface (see Fig. 5.9(e)) shows that: (1) The width of the interface is large, spanning about 
10 — 15 atomic spacings and the interface is wet by a buckling phase made up of triangular 
groups of nine- layered solid shifted in a direction normal to the walls. (2) The interface between 
rii layered crystal and rii — \ layered smectic contains a dislocation -^with Burger's vector in the 
y- direction which makes up for the difference in the number of layers. Each band of width s is 
therefore held in place by a dislocation-anti-dislocation pair (Fig. 5.9). In analogy with classical 
nucleation theory[148, 149], the free energy F;, of a single band can be written as 



A comparison of Fig. 5.9 (c) with Fig. 1 of Rcf.[117] leads us to speculate whether a two-dimensional 
version of the "buckled" phase may in-fact wet the solid-smectic interface thereby reducing its energy. 



e. 



= -6Fs + Ec + ^b^K^ log — 




(5.5) 
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Figure 5.12: Dislocation probabilities of a 65 x 10 system are plotted as a function of strain 
starting from a solid commensurate with the interwall spacing Ly. The + symbols denote dislocation 
probabilities for only those Burger's vectors which have components perpendicular to the walls. The 
corresponding bond-breaking moves are of the type depicted by the uppermost plot in the inset panel. 
Arrows show the directions of the bond-breaking moves. On the other hand * symbols denote other 
dislocation probabilities corresponding to the other two types of moves shown in the inset. 

where b — ay/2 is the Burger's vector, SF the free energy difference between the crystal and the 
smectic per unit length and Ec the core energy for a dislocation pair. Bands form when dislocation 
pairs separated by s > ^y^K^/SF arise due to random fluctuations. To produce a dislocation 
pair a large energy barrier of core energy Ec has to be overcome. Though even for very small 
strains Ed the free energy becomes unstable the random fluctuations can not overcome this 
large energy barrier within finite time scales thereby suppressing the production of — 1 layered 
smectic bands up to the point of e^. In principle, if one could wait for truely infinite times the 
fluctuations can produce such dislocation pairs for any non-zero Ed though the probability for such 
productions exp(— /3£'c) si's indeed very low. Using a procedure similar to that used in Ref.[71], 
we have monitored the dislocation probability as a function of 77 (Fig. 5. 12). Not surprisingly, the 
probability of obtaining dislocation pairs with the relevant Burger's vector increases dramatically 
as 7/ — r/ci and artificially removing configurations with such dislocations suppresses the transition 
completely. Band coalescence occurs by diffusion aided dislocation "climb" which at high density 
implies slow kinetics. Due to this diffusive nature the size of smectic band L{t) scales as L ~ \/t. 
Throughout the two-phase region, the crystal is in compression and the smectic in tension along 
the y direction so that a is completely determined by the amount of the co-existing phases; 
orientation relationships between the two phases being preserved throughout. Again the amount 
of solid or smectic in the system is entirely governed by the strain value e. This means that for 
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a given strain tine amount of solid and smectic and therefore the amount of stress a is entirely 
determined by the value of strain. This explains the reversible[150] plastic deformation in Fig. 



5.4 Mean Field Results: The Reversible Failure Transition 

The failure of a commensurate solid under tensile strain imposed in the manner discussed in the 
previous section, comes about through the nucleation of smectic bands within the solid. Monte- 
Carlo simulations show, at half-integral x where the local internal strain Sd becomes discontinuous, 
p(r) at nearest neighbour sites overlap along the x-direction, parallel to the walls, generating 
smectic bands. The stress associated with Sd vanishes at these points and the solid fails under 
tension. In this section we shall show, using simple density functional[148, 151] arguments, that 
the phase transition and the consequent tensile failure (a smectic cannot support stress parallel 
to the smectic layers) is brought about by this overlap in the local density. Since mechanical 
failure in our system is a consequence of a phase transition, it is reversible — as the strain is 
reduced back to zero, the stress also vanishes and the perfect triangular lattice is recovered[152]. 

Within density functional theory[151], the excess grand potential of a non-uniform liquid con- 
taining a density modulation p(r) over the uniform liquid of density pi is given by, 



Here 5p{v) — p(r) — pi and C(r) is the direct correlation function of the uniform liquid[109]. 
A functional minimization of the free energy yeilds the following self-consistency equation for the 
density: 



In principle one should solve the above equation within the constraints imposed by the walls 
and obtain the equilibrium p(r). Substitution of this p(r) into Eqn.5.6 gives the equilibrium free 
energy and phase transitions. While we intend to carry out this procedure in the future, we must 
point out that for the present problem, this is complicated by surface terms and anisotropic, 
external, fields which are difficult to incorporate. In this chapter we shall take a much simpler 
route in exploring the various conditions for the solid -smectic transition given the nature of the 
PGi (the order parameters) obtained from our simulations. 



5.8. 




(5.6) 




(5.7) 
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One may expand, therefore, the logarithm of the local density profile logp(r) in a Fourier 
series[148] around a lattice point at the origin, to get, 

p{r)=Xexp ^2Co5^pGiCos(Gi.r)j (5.8) 

where Cq is a constant, of order unity, denoting the Fourier transform of the direct correlation 
function calculated at a q-vector corresponding to the smallest RLV set of the solid. We have 
kept contributions only from this set. 

For a perfect triangular lattice, the RLV's are Gi = y^, G2 = £^cos(|) +?j^sin(|) and 
G3 = x^cos{^) — y^sm(^), where dy = ^a^. Using these relations and the fact that in 
the presence of confining walls the Fourier amplitudes denoting solid order are virtually constant 
upto the transition and = Pca 7^ PGi. Eqn.5.8 gives [153], 

p(r) =A/'exp{Co(2pG, +4pGj}exp ^-^Cq {(2pG, + PgJi/' + SpG.a;'} j (5.9) 

Clearly the density profile is Gaussian, of the form, p(r) ~ exp{—y'^ /2ay — x'^/2a1). Therefore, 
the spreads of density profile in x and |/-directions are given by and ay respectively, with 

in the absence of walls, pci = PG2 making — ay, i.e. the density profile comes out to be 
symmetric in both directions, as expected for the bulk triangular solid. The presence of walls make 
ay < a^ making the density profile elliptical with larger spread in x-direction, the direction parallel 
to the walls. Two neighbouring density profiles will overlap to form a smectic if ax > a^. This 
leads us to the definition of a measure of overlap Oi — {ax/oLx) ■ The condition > 1 is then the 
Lindemann criterion for nucleation of the smectic phase. Remembering — ao{ni — l)/(x — 1), 
we get, 

1 1 y - 1 
Oi = ^ — -. (5.12) 

Whenever, pc^ i.e. with the loss of solid order Oi diverges although ay remains finite, 
since pGi 7^ in presence of the walls. This indicates a solid-smectic transition. However, even 
before ~^ the quantity A = and therefore Oi shows large jumps at those internal 
strain [sd) values where x becomes half-integer. It is interesting to note that, at these points 
Ed has discontinuities and the system fails[152]. This shows that the mode of failure predicted 
by our theory is through a solid-smectic transition. The fact that pgi remains non-zero even at 
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Figure 5.13: For a 10-layered solid with Ly commensurate with the initial strainless triangular 
structure at 77 = .85 overlap term A is plotted as a function of external strain e. Density profile overlap 
shows a jump increase at strains e > .1, the failure strain value[152]. 

very small densities, due to the confinement from the walls, gives rise to the density modulations 
in the confined liquid. 

We have shown therefore that the overlap in the density profiles may be used as an "order 
parameter" for the solid to smectic transition. We show below that jumps in this order parameter 
tantamounts to mechanical failure of the solid. 

We begin with studying the overlap A as a function of external strain e. For specificity, 
we start from a triangular solid of packing fraction -q = .85 with Ly commensurate with a 
rii = 10 layered solid. With increasing strain initially the overlap A reduces due to increased 
separation {a^) between neighbouring lattice points. But above a strain (e) of about 10%, x 
reaches the half-integral mark and A shows a discontinuous increase, indicating large overlap 
between neighbouring density profiles along the wall ; indicating a solid to smectic transition 
(Fig. 5. 13) at the failure strain e*. With further increase in strain the overlap reduces, again due 
to increased separation between neighbouring lattice points. At higher strains the smectic melts 
into a modulated liquid due to increased fluctuations connected with the reduced density[152]. 

We have performed this calculation for various Ly values commensurate with starting triangular 
solids of = 2, 3 . . . 20 layers at packing fraction i] = .85. We found out the failure strains e* 
at each Ly and plotted them in Fig. 5. 14 as a function of Ly. This clearly shows that the failure 
strain reduces with increase in Ly. This demonstrates the fact, derived earlier from Monte-Carlo 
simulations[152], that thinner (smaller Ly) strips are stronger! 

In Fig. 5. 15 we have plotted the overlap term A with increasing interwall separation Ly at 
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Figure 5.14: Failure strains e* for various interwall separations Ly confining = 2 — ^ 20 layered 
triangular strips at 77 = .85 is plotted. Failure strain decreases with increase in Ly. 

r] = .85. The jumps, as usual, indicate failure strains corresponding to discontinuities in the 
internal strain Sd at half-integral values of x- The plot shows that the amount of overlap at 
the failure strains e* reduces with increasing Ly indicating that at large interwall separations the 
system starts to behave as a bulk solid and more conventional modes of failure viz. through 
formation and interaction of cracks and twin boundaries starts becoming active. 



5.5 Mean Field Results: The Equilibrium Phase Diagram 

We now calculate the equilibrium phase diagram of the confined two dimensional system in i] — Ly 
plane. For this purpose we utilise the usual common tangent construction method to extract the 
phase diagram in the canonical ensemble from the free energies of the solid and fluid, the only 
two unambiguous stable phases of hard disk system in thermodynamic limit. Our computation 
will, evidently, ignore the other possible phases like smectic and modulating liquid. Physically, 
confinement always induces density modulations and therefore in phase diagram we understand 
that the liquid is always a modulated liquid for small channel widths. The phenomenological 
equation of state due to Santos et. a/.[154] agrees with bulk fluid simulations of hard disks. The 
corresponding excess free energy per particle is 

(2,. - 1) log (1 - (2v, -1)^)- log(l - 

?7c being the packing fraction of hard disk solid at close-packed limit. The total free energy of 
the fluid will have, moreover, a contribution of free particle free energy per particle (logp— 1). 
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Figure 5.15: Overlap A is plotted for a system at rj = .85 with changing interwall separation L 
Amount of smectic overlaps A at failures reduces with increasing Ly. 



Then the fluid free energy per unit volume becomes Tpi — p{J^San + logp — 1)- To incorporate 
the effect of structureless walls on the fluid we have incorporated the surface tension of a hard 
disk fluid obtained from scaled particle theory (SPT) [137], 

1.1 = -TTTT^pd + UPd (5.14) 



2(1-77)2'^ 2' 



where the pressure is given by [154] 



2 



= I 1-277+ (^j (2r;c-l)) (5.15) 

Then the liquid free energy per unit volume becomes jFpj = p{J-'san + logp — 1) + Ipi/Ly. 

We calculate the solid free energy of hard disk system from free volume theory. The surface 
tension of solid in presence of walls is given by[137] 



1 d 

y/Saia-d) [ 2 



Hence the solid free energy per unit volume is J-'cr = ^t{v^ x)~^lcrl -^Tiv^ x) is the FNFVT 
free energy as described in Sec. 5. 3. For every wall to wall separation Ly we vary the packing 
fraction i] and pick up the smaller free energy at each r] to obtain F = m\f\[J^cr , ^fi) in Fig. 5. 16. 

Then we use common tangent construction over this free energy F to obtain densities of 
coexisting phases at equal chemical potential. Within this theory, at densities lower than the 
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Figure 5.16: Plot of free energy used for constructing the phase diagrann. The dotted curve, dashed 
curve and solid curve denote Tpu ^Cr and F respectively. 



Figure 5.17: Phase diagram. White region denotes solid phase. The darkest region at high densities 
denote regions of phase space inaccessible to the system due to overlapping hard disks. The darker (aqua 

green) region at low densities denote fluid phase. All other regions are two phase coexistence regions. 
The bold (blue) dashed lines denote x = ra; and the faint (green) dashed lines denote half integral x- 

coexisting fluid density the thermodynamic phase is stable fluid and at densities higher than the 
coexisting solid density it is stable solid. At intermediate densities fluid and solid coexist. The 
corresponding phase diagram is given in Fig. 5. 17. 

From the phase diagram it is evident that solid phase at lower density gets more and more 
instable as we increase the wall to wall separation. If we start with solids commensurate with 
Ly, i.e. X = n/ at packing fraction rj = .85 and start reducing r] keeping Ly fixed, we can find 
out from the phase diagram the deviatoric strains at which one hits the two-phase boundary as 
well as the strains where one reaches at half-integral values of %. These have been plotted in 
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Fig.5.18. 



0.14 



0.12 




0.1 










0.06 




0.04 - 




0.02 L 



7 



8 



9 



TO rr 



12 13 14 15 




Figure 5.18: The data points denoted by * show the deviatoric strain values at which x = rt; — 1/2 
whereas the symbols + denote the arrival of two phase coexistance, i.e. failure predicted by this theory. 
Lines are guide to eye. Data have been extracted from Fig.5.17. 

Both the quantities show a decrease in failure strain £^ with increase in width of the strip 
Ly supporting the simulation result, narrower strips are stronger. Our simple theory therefore 
predicts a first order solid-fluid transition as a function of Sa and failure of the solid as the density 
enters the region of solid-fluid coexistence at larger critical strain £^ for smaller wall to wall 
separations. However, details like the density modulation, effects of asymmetry in density profile, 
vanishing displacement modes at the walls and most importantly nucleation and dynamics of 
misfit dislocations crucial to generate the smectic band mediated failure observed in simulations 
are beyond the scope of this theory. Also, the effect of kinetic constraints which stabilize the solid 
phase well inside the equilibrium two phase coexistence region is not captured in this approach. 
We believe, nevertheless that this equilibrium calculation may be used as a basis for computations 
of reaction rates for addressing dynamical questions in the future. 

5.6 Fluctuations and Destruction of Order 

One of the key definitions of a solid states that a solid, as opposed to a liquid, is a substance 
which can retain its shape due to its nonzero shear modulus. Going by this definition, a QlD 
solid confined within planar, structureless walls is not a solid despite its rather striking triangular 
crystalline order as well as an apparently solid- like structure factor. Indeed, the shear modulus of 
the confined solid at = 0.85 is zero, though the corresponding system with periodic boundary 
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Figure 5.19: Shear stress vs. shear strain at = 0.85. A system of 43x50 2D hard disks simulated 
with periodic boundary conditions gives a shear modulus fj. = 393.42 with an error within 1%. On the 
other hand when a 40 x 10 triangular lattice of hard disks is confined within a commensurate channel, 
that fits 10-layers of lattice planes, the shear modulus drops drastically to /i = 0! 

condition show large shear modulus (See Fig. 5. 19). This is a curious result and is generally true 
for all values of 4 < < 25 and r] investigated by us. 

To understand the nature and amount of fluctuations in the confined QlD system we calculate 
correlation between displacement fields along the channel, < (u^{x) — u^(0))'^ > for a layer of 
particles near a boundary. The nature of the equilibrium displacement correlations ultimately 
determines the decay of the peak amplitudes of the structure factor and the value of the equi- 
librium elastic moduli [148]. In one dimension < (u^{x) — u^{0)Y >~ x and in two dimensions 
< {u^{x) — u^{0)y >~ ln(x). In the QlD system it is expected that for small distances, dis- 
placement fluctuations will grow logarithmically with distance which will crossover to a linear 
growth at large distances [155, 156]. We calculate this quantity averaged over 10, 30, 50, 100 
configurations equlibrated over 10^ Monte- Carlo (MC) steps and separated by lO'^ MC steps. 
We compare the results obtained from a 5000 x 10 hard diks solid at r/ = 0.75 with periodic 
boundary conditions (PBCs) and with hard channel confinement (Fig. 5. 20). The structure factor 
for this system is apparently solid like with prominent triangular order. The shear modulus is 
vanishingly small for the confined system and non- zero for a system with PBCs. We calculate 
fluctuations averaged over 10, 30, 50, 100 configurations equlibrated over 10^ Monte- Carlo 
(MC) steps and separated by 10'^ MC steps. The fluctuation of displacement field for the system 
with PBC reduces and converges with the increase in number of configurations over which aver- 
agings are done. However, for the confined solid displacement fluctuations continue to increase 
with number of configurations over which the averaging is done. This clearly shows that as soon 
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as a solid gets confined, even if the confinement length scale is commensurate with that set by 
the system density, the solid starts to behave like a liquid with zero shear modulus and linearly 
increasing displacement correlations. Density modulations introduced by the hard walls seems to 
destabilize long ranged order in the small channel. This is reminiscent of large potential strength 
limit of laser induced transitions [89, 94, 157] as discussed in the last chapter. However, when 
we compare the results obtained from a 5000 x 10 hard disks solid at 77 = 0.85 with PBCs and 
hard channel confinement (Fig. 5. 21) we obtain a different result. Fig. 5. 21 clearly shows that with 
increase in number of configurations over which the averagings are done, fluctuation reduce and 
converge to a small number. This is consistent with the solid like structure factor but inconsistent 
with the vanishing shear modulus of confined system. This result is particularly surprising since 
it can be shown[155] that all QID solids whether with PBCs or not and for any density 77 would 
have displacement correlations which increase linearly with system length. The exact value of 
this crossover length has nonuniversal prefactors which depend on boundary conditions, density, 
nature of interactions etc. If the crossover length turns out to be larger than the system size, 
then this increase will not be observed and the properties of the system will be anomalous and 
inconsistent. The fact that displacement fluctuations saturate at high densities as seen in our 
simulations suggests that this may, in fact, be the reason behind the puzzling behaviour. An 
accurate calculation of the crossover length is required to settle this question conclusively - a 
task which is non-trivial due to the nonuniversal nature of this number. 

A separate calculation of Young's moduli (response to elongational strain) on a commensurate 
40 X 10 hard disk confined solid at 77 = 0.85 shows that = 1361 and Yy = 1503 within 3% 
error. Young modulus in the longitudinal direction is smaller than that in the direction transverse 
to the confinement and both these values are larger than the Young modulus of the system under 
PBC {Y — 1350). This implies that the non-hydrodynamic component of the stress a^x — cfyy 
is non-zero for non zero strains as shown in Fig.5.8. Therefore even if we choose to regard this 
QID solid as a liquid, it is quite anomalous since it trivially violates Pascal's law which states 
that the stress tensor in a liquid is always proportional to the identity! Lastly, commensurability 
seems to affect strongly the nature and magnitude of the displacement fluctuations which increase 
dramatically as the system is made incommensurate. Similar behaviour has been recently noticed 
by A. Ricci et. al. [155] for a confined soft disk system. They also noticed that the ID structure 
factor of their system showed liquid like behaviour. Here we note that in the work done by A. 
Ricci et. al. [155], they used a soft core interaction potential (1/r^^) and a soft system- wall 
potential (l/r^°) and studied the shear modulus and structure factor at density p = 1.05, very 
close to bulk phase transition at ksT = 1. In that system it is difficult to determine whether the 
system is commensurate or not, since the long ranged soft wall potential would tend to distort 
the triangular lattice so that the equilibrium state would not be a perfect triangular solid. Our 
system has the advantage of having much simpler ground states which makes these distinctions 
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unambiguous. 

5.7 Conclusion 

In this chapter we explored some of the strange and anomalous properties of QID solids which 
are only a few atomic dimensions wide in one direction. We have found that these properties 
persist even if the width is increased by a large amount and the approach to bulk behaviour is 
slow. The bulk limit is also approached in an oscillating manner with commensurability playing a 
very important role. 

What impact, if any, do our results have for realistic systems? Apart from constrained hard 
sphere colloids[96] where our results are directly testable, a similar fracture mechanism may be 
observable in experiments on the deformation of mono-layer nano beams or strips of real materials 
provided the confining channel is made of a material which is harder and has a much smaller 
atomic size than that of the strip[113, 114]. The effect of elasticity and corrugations of the walls 
on the fracture process, as well as it's dynamics, are interesting directions of future study. The 
destruction of long ranged solid like order should be observable in nano wires and tubes and may 
lead to fluctuations in transport quantities [158]. 

In the next chapter we shall study the transport of heat accross a QID solid and try to find out 
the effect of the reversible failure transition in the heat transport coefficient. 
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Figure 5.20: < {u^{x) — u^(0))^ > /cP fluctuations within a single lattice plane (line) is averaged 
over equilibrated configurations of 5000 x 10 systenn of hard disks at r] = 0.75 and simulation box 
dimensions which are commensurate with the lattie parameter of the bulk (unconfined) system. The 
upper panel shows the case for a system with periodic boundary condition while the lower panel shows 
the same for QID confinement using hard walls. For the confined system, fluctuations increase with 
increase in the number of configurations over which the displacement correlations are averaged without 
showing any sign of convergence. Fluctuations in the system with periodic boundary conditions, on the 
other hand, appear to saturate. 
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Figure 5.21: The same quantities as in Fig. 5. 20 averaged over equilibrated configurations of 5000 x 
10 system of hard disks at = 0.85 and a simulation box size that is commensurate with the lattie 
parameter of the bulk (unconfined) system. With increase in the number of configurations over which 
averaging is done, fluctuations reduce and converge, both for a system with periodic boundary condition 
and for the confined system in channel. 
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6 Heat Conduction in Confined Solid 



Jose Arcadio Buendia ventured a murmur: "It's the largest diamond in the 
world". "No", the gypsy countered, "It's ice". - G. G. Mdrquez 



In the previous chapter [152] we have observed that the properties of a solid that is confined 
in a narrow channel can change drastically for small changes in applied external strain. This 
was related to structural changes at the microscopic level such as a change in the number of 
layers of atoms in the confining direction. These effects occur basically as a result of the small 
(few atomic layers in one direction) dimensions of the system considered and confinement along 
some direction. A similar layering transition, in which the number of smectic layers in a confined 
liquid changes discretely as the wall-to-wall separation is increased, was noted in [132, 133]. 
Both [132, 152] look at equilibrium properties while [133] looks at changes in the dynamical 
properties. An interesting question is as to how transport properties, such as electrical and 
thermal conductivity, get affected for these nanoscale systems under strain. These questions are 
also important to address in view of the current interest in the properties of nanosystems both 
from the point of view of fundamentals and applications [20, 159, 160]. 

In this chapter we consider the effect of strain on the heat current across a two-dimensional 
"solid" formed by a few layers of interacting atoms which are confined in a long narrow channel. 
We note here that, in the thermodynamic limit it is expected that there can be no true solid 
phase in this quasi-one-dimensional system. However for a long but finite channel, which is our 
interest here, and at high packing fraction the fluctuations are small and the system behaves like 
a solid. We will thus use the word "solid" in this sense. 

In previous chapter [152] the anomalous failure, under strain, of a narrow strip of a two 
dimensional solid formed by hard disks confined within hard walls [ see Fig. 6.1 ] was studied. 
Sharp jumps in the stress-strain were observed. These were related to structural changes in the 
system which underwent transitions from solid-to-smectic-to-modulated liquid phases [152, 153]. 
In the present chapter we study changes in the thermal conductance of this system as it undergoes 
elastic deformation and failure through a layering transition caused by external elongational strains 
applied in different directions. 

The calculation of heat conductivity in a many body system is a difficult problem. The Kubo 
formula and Boltzmann kinetic theory provide formal expressions for the thermal conductivity. 
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Figure 6.1: A solid with a triangular lattice structure formed by hard disks confined between two 
structureless walls at y = and y = Ly. The walls are maintained at two different temperatures. The 
lattice parameters of the unstrained solid are denoted by a% and a^. Elongational strains can be imposed 
by rescaling distances either in the x or y directions and the lattice parameters change to and ay. 



In practice these are usually difficult to evaluate without making drastic approximations. More 
importantly a large number of recent studies [161-164] indicate that the heat conductivity of 
low-dimensional systems infact diverge. It is then more sensible to calculate directly the heat 
current or the conductance of the system rather than the heat conductivity. In this chapter 
we propose a simple-minded calculation of the heat current which can be expected to be good 
for a hard disk (or hard spheres in the three dimensional case) system in the solid phase. This 
reproduces some qualitative features of the simulations and gives values for the current which are 
of the correct order of magnitude. 

The organization of this chapter is as follows. In Sec. (6.1) we explain the model and present 
the results from simulations. In Sec. (6.2) we derive a simple formula for heat current in a hard- 
sphere system and evaluate it approximately. We conclude with some discussions in Sec. (6.3). 

6.1 Results from Simulations 

We consider a two dimensional system of hard disks of diameter d and mass m which interact with 
each other through elastic collisions. The particles are confined within a narrow hard structureless 
channel [see Fig. 6.1]. The hard walls of the channel are located at y = and y = Ly and we 
take periodic boundary conditions in the a:— direction. The length of the channel along the 
x— direction is L^. and the area is ^ = x Ly. The confining walls are maintained at two 
different temperatures { T2 at y = and Ti at y = Ly ) so that the temperature difference 
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AT — T2 — Ti gives rise to a heat current in the y-direction. Initially we start with channel 



a triangular lattice structure. We then study the heat current in this system when it is strained 
(a) along the x— direction and (b) along the |/— direction. 

We perform an event-driven collision time dynamics [165] simulation of the hard disk system. 
The upper and lower walls are maintained at temperatures Ti = 1 and T2 — 2 (in arbitrary 
units) respectively by imposing Maxwell boundary condition [161] at the two confining walls. 
This means that whenever a hard disk collides with either the lower or the upper wall it gets 
reflected back into the system with a velocity chosen from the distribution 



where Tw is the temperature (Ti or T2) of the wall on which the collision occurs. During 
each collision energy is exchanged between the system and the bath. Thus in our molecular 
dynamics simulation, the average heat current flowing through the system can be found easily 
by computing the net heat loss from the system to the two baths (say Qi and Q2 respectively) 
during a large time interval r. The steady state heat current from lower to upper bath is given 
by (/) = limt-.ooQi/'T — — limt_»oo <52/t. In the steady state the heat current (the heat flux 
density integrated over x) is independent of y. This is a requirement coming from current 
conservation. However if the system has inhomogeneities then the flux density itself can have 
a spatial dependence and in general we can have j — j{x,y). In our simulations we have also 
looked at j{x, 0) and j{x, Ly). 

Note that the relevant scales in the problem are: ksT for energy, d for length and = 
■\/m(P /kBT for time. We start from a solid commensurate with its wall to wall separation and 
follow two different straining protocols. In case (a) we strain the solid by rescaling the length in 
the x— direction and the imposed external strain is e,j,j, = (L^ — L°.)/L°. In case (b) we rescale 
the length along the |/— direction and the imposed strain is tyy = {Ly — L°)/L°. 

The only thermodynamically relevant variable for a hard disk system is the packing fraction 
7] = TiNcP /AA. For a close packed solid with periodic boundary condition this value is about r]c = 
0.9069. On the other hand for a confined solid having Ny number of layers rjc = 7rNy/{2\/3(Ny — 
1) + 4) and for a 10- layered solid r]c = 0.893. In our simulations we consider initial values of 
7] for the solid to be close to 77^. The channel is "mesoscopic" in the sense that it has a small 
width with Ny = 10 layers of disks in the y— direction (in the initially unstrained solid). In the 
direction the system can be big and we consider A^^; = 20, 40, 80, 160 number of disks in 
the x— direction. In collision time dynamics we perform 10"^^ collisions per particle to reach steady 
state and collect data over another 10^ collisions per particle. All the currents calculated in this 
study are accurate within error bars which are less than 3% of averaged current. 



dimensions L\ 



an 



d Ly such that the system is in a phase corresponding to a unstrained solid with 
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Figure 6.2: Plots obtained by superposition of 500 steady state configurations of a portion of 
40 X 10 system taken at equal time intervals. Starting from rj = 0.85 imposition of strains (a) exx = 0.1, 
(b) exx = 0.15 gives rise to these structures. The colors code local density of points from red/dark 
(high) to blue/light (low). In (a) one can see a 9-layered structure nucleated within a 10-layered solid. 
The corresponding structure factor identifies this to be a smectic [152]. In (b) the whole system has 
transformed into a 9-layered smectic. 



Let us briefly recapitulate some of the equilibrium results for stress-strain behavior obtained in 
the last chapter. As the strain exx is imposed, the perfectly triangular solid shows rectangular 
distortion along with a linear response in strain versus stress behavior. Above a critical strain 
(^exx ~ 0.1) one finds that smectic bands having a lesser number of layer nucleate within the 
solid [this can also be seen in Fig. (2a), which is for a nonequilibrium simulation]. This smectic 
is liquid-like in a;— direction (parallel to the walls) and has solid-like density modulation order in 
I/— direction (perpendicular to the walls). With further increase in strain the size of the smectic 
region increases and ultimately the whole system goes over to the smectic phase at exx ~ 0.15 
[Fig. (2b)]. At even higher strains the smectic melts to a modulated liquid [152, 153]. The 
modulated liquid shows typical liquid like ring pattern corresponding to average inter- particle 
separation above the smectic like ID density modulation peaks in the structure factor [152, 166]. 
This layering transition is an effect of finite size in the confining direction. Similar phase behaviors 
have been observed in experiments on steel balls confined in quasi ID [134]. We note that to fit 
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Figure 6.3: Plot of jy (in units of kBT/Tgd) versus exx foi' different lengths of the channel. The 
width of the channel is Ny = 10 layers. Starting packing fraction is r/ = 0.85. The solid line shows the 
theoretical prediction of dependence of the heat current on strain [see Sec. (6.2)]. 

a Ny layered triangular solid within a channel of width Ly we require 

Ly = ^al{Ny-l) + d. (6.2) 
This enables us to define a fictitious number of layers 

of triangular solid that can span the channel where a is the lattice parameter at any given 
density. The actual number of layers that are present in the strained solid is Ny = I{x) where 
the function I{x) gives the integer part of x- For confined solids the free energy has minima at 
integer values of x and maxima at half-integral values [152, 153]. The difference in free-energy 
between succesive maxima and minima gradually decreases with increasing Ly. Thereby the 
layering transition washes out for ni = 25 layered unstrained solid[152]. Up to this number of 
layers, a triangular solid strip confined between two planar walls fails at a critical deviatoric strain 
ejj ~ '^/Ny, i.e. smaller strips fail at a larger deviatoric strain (e^ = e^x — ^yy)- 

We now present the heat conduction simulation results for the two cases of straining in x and 
y directions. 

(a) Strain in x direction - In Fig. 6.3 we plot the heat current density jy calculated at different 
values of the strain txx- Starting from the triangular lattice configuration, we find that the 
heat current decreases linearly with increase in strain. At about the critical strain txx ~ 0.1 we 
find that the heat current begins to fall at a faster rate. This is easy to understand physically. 
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Figure 6.4: is the local number of layers averaged over 10'^ steady state configurations for a 

system of 40 x 10 hard disks. A starting triangular lattice of r/ = 0.85 is strained to exx = 0.118 and 
the data collected after steady state set in. Also shown is the local heat current jy{x). The regions 
having lower number of layers conduct less effectively. 

At the onset of critical strain smectic bands, which have lesser number of particle layers, start 
nucleating (Fig. 6.2). These regions are much less effective in transmitting heat than the solid 
phase and the heat current falls rapidly as the size of the smectic bands grow. At about the 
strain value exx ~ 0.15 the whole system is spanned by the smectic. Beyond this strain there is 
no appreciable change in the heat current. The solid line in Fig. 6.3 is an estimate from a simple 
analysis explained in Sec. (6.2). 

In Fig. 6.4 we plot the local steady state heat current jy{x) for a system of 40 x 10 particles at 
a strain exx = 0.118 i.e. at a strain corresponding to the solid-smectic phase coexistence. At this 
same strain the number of layers averaged over 10^ configurations have been plotted. It clearly 
shows that the local heat current is smaller in regions with smaller number of layers. This is the 
reason behind getting a sharp drop in average heat current after the onset of phase coexistence. 

(b)Strain in y direction - Next we consider the case where, again starting from the density 
r] = 0.85, we impose a strain along the y— direction. As shown in Fig. 6.5, the heat current jy 
now has a completely different nature. The initial fall is much steeper and has a form different 
from the linear drop in Fig. 6.3. The approximate analytic curve is explained in Sec. (6.2). At 
about eyy ~ 0.1 we see a sharp and presumably discontinuous jump in the current. At this 
point the system goes over to a buckled phase (Fig. 6.6b) in which different parts of solid (along 
X- direction) are shifted in y-direction by small displacement to cover extra space between the 
walls[118-120]. A further small strain induces a layering transition and the system breaks into 
Ny = 11 layered solid and Ny = 10 layered highly fluctuating smectic like regions. At even higher 
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Figure 6.5: Plot of jj, (in units of kBT/Tgd) versus eyy for different channel lengths. The channel 
width is Ny = 10 layers. The starting packing fraction is 77 = 0.85. The jump in current occurs at 
the strain value where the number of layers in the y— direction increases by one and the system goes to 
a smectic phase. The solid line shows the theoretical prediction of dependence of the heat current on 
strain [ see Sec. (6.2) ]. 



strains {tyy ~ 0.2) the whole system eventually melts to a A^^ = 11 layered smectic phase. The 
phase behavior of this system is interesting and will be discussed in detail elsewhere[166]. Unlike in 
the case with applied strain in the x— direction, in the present case the buckling-layering transition 
is very sharp. Even though the overall density has decreased, due to buckling and increase in 
number of layers in the conducting direction, there is an increase in the energy transferring 
collisions and hence the heat current. The plots in Fig. 6.6 show the structural changes that 
occur in the system as one goes through the transition. 

We find in general that the heat current along any direction within the solid follows the same 
qualitative features as the stress component along the same direction. This can be seen in 
Fig. 6. 7 where we have plotted jy versus exx for two starting densities of solids 77 = 0.85, 0.89. 
In the inset we show the corresponding —o-yy versus exx curves and see that they follow the same 
qualitative behavior as the heat current curves. The reason for this is that microscopically they 
both originate from interparticle collisions. Infact the microscopic expressions for the total heat 
current [ see Eq. (6.11) in Sec. (6.2) ] is very similar to that for the stress tensor component, 
with an extra velocity factor. The stress tensor is given by: 

Aa., = -5;(»„>f) + ^/M!ii)M\ , (6.3) 

i i<j \ ^^'^ I 
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Figure 6.6: Plots obtained by superposition of 500 steady state configurations of a portion of 40 x 10 
system taken at equal time intervals. Starting from rj = 0.85 imposition of strains (a) Eyy = 0.05, (b) 
eyy = 0.1, (c) eyy = 0.12 gives rise to these structures. The colors code local density of points from 
red/dark (high) to blue/light (low), (a) Solid phase, (b) A mixture of 10-layered solid and a buckling 
phase, (c) An 11-layered solid in contact with 10-layered smectic like region. 

where {x'^,u°'} refer to the a-th component of position and velocity of the i^^ particle, rfj = 
Sa(^ij)^ and (pivij) is the interparticle potential. For a hard disk system, ^^^^^ can be replaced 
by —ksTdirij — d). Also in equilibrium we have {mufu^} = ksTdap and hence the stress tensor 
becomes: 




Using collision time simulation it is easier to evaluate the stress tensor in the following way. We 
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Figure 6.7: Plot of jy (in units of ksT/Tsd ) versus €xx foi" two different stsrting vslues of the 
packing fraction. corresponds to a starting value of rj = 0.89 while + is for rj = 0.85. In both 
the cases the initial solid size was 80 x 10. The inset shows corresponding plots of —(Tyy (in units of 
kBT/d?) versus e^x- Notice that stress-strain curve has the same qualitative profile as the jy versus Exx 
curve. 

can rewrite Eq. (6.3) as 

i<j 

We use the idea that (. . .) can be replaced by a time average so that from Eq. (6.3) we have 

Now note that during a collision we have / dtfj^j = Ap^j where Apij is the change in momentum 
of i^^ particle due to collision with j^^ particle. It can be shown that Apij = —(uij.fij)fij where 
fij = Tij/vij and Uij = Ui — Uj and r^, Ui are evaluated just before a collision. This change in 
momentum occurs for a single pair of particle during one collision event. To get the stress tensor 
we sum over all the collision events in the time interval r between all pairs of particles. Therefore 
in collision time dynamics we get the following expression for the stress tensor, 

Aa^p = -NkBTS^p + Jim ^ J] E ^P^4 ' (6-4) 

Tc i<j 

where V denotes a summation over all collisions in time r. 
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6.2 Analysis of Qualitative Features 

We briefly outline a derivation of the expression for the heat flux. For the special case of a hard 
disk system this simplifies somewhat. We will show that starting from this expression and making 
rather simple minded approximations we can explain some of the observed results for heat flux 
as a function of imposed external strain. 
We consider a system with a general Hamiltonian given by: 

^=E[^ + ^(rO] + ^E'^(^^^)' (6-5) 

where V"(rj) is an onsite potential which also includes the wall. To define the heat current density 
we need to write a continuity equation of the form: de{r,t)/dt + dja{r,t)/dxa = 0. The local 
energy density is given by: 

e{r,t) — (5(r — rj)/ij where 



Taking a derivative with respect to time gives 

I = -£-J2^{r-v,)h,ut + J25iv-v,)h, (6.6) 

i i 

where = '^iS{r — rj)/ijUj is the convective part of the energy current. We will now try to 
write the remaining part given by as a divergence term. We have 



i 

i * 



2 



where f^j — —d(f){rij)/dxf. Using the equation of motion muf — —dV/dxf + ^j^^f^ we get 

= jE^(r-r')(/x-/x) • (6-8) 



2 . 



With the identification W'^ = —dj^/dx°' and using fjj = — fjj we finally get: 

fM = ^ E ^(^" - n '^(^'^ - + (6-9) 



2 
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where 6{x) is the Heaviside step function. This formula has a simple physical interpretation. First 
note that we need to sum over only those i for which xf > x". Then the formula basically gives 
us the net rate at which work is done by particles on the left of x°' on the particles on the right 
which is thus the rate at which energy flows from left to right. The other part, , gives the 
energy flow as a result of physical motion of particles across x°'. Let us look at the total current 
in the system. Integrating the current density over all space we get: 



2 . .,. 

= (6.10) 

Including the convective part and taking an average over the steady state we finally get: 

i 

We note that for a general phase space variable A{{xi^Ui^) the average (A) is the time average 
lim,^oo(l/r) g dtA{{xi{t),Ui{t)\). 

Finding the energy current for hcird disk system: The energy current expression 
involves the velocities of the colliding particles which change during a collision so we have to be 
careful. We use the following expression for (/^): 



i,3+i 

- i-^r4E4(/^-f-/>') (6-12) 

Now if we integrate across a collision we see that J dt{fij.Ui) gives the change in kinetic energy 
of the i^^ particle during the collision while J dt{fji.Uj) gives the change in kinetic energy of the 
j^^ particle. Hence we get 
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where we have used the fact that for elastic collisions Ai^j = — AXj and denotes a sum- 
mation over all collisions, in the time interval r, between pairs {ij}- The time interval between 
successive collisions between i*-^ and j^^ particles is denoted by and the average ( ... )c in 
the last line denotes a collisional average. Thus (rjj)c = liniT-^oo T/-^jj(T), where Nij{T) is the 
number of collisions between i^^ and j*^ particles in time r. For hard spheres the convective 
part of the current involves only the kinetic energy and is given by ( ) = ( {muf/2)uf ). 
Using these expressions we now try to obtain estimates of the heat current and its dependence 
on strain in the close packed limit where the system looks like a solid with the structure of a 
strained triangular lattice. 

In the close packed limit the convection current can be neglected and we focus only on the 
conductive part given by (/'^) = {I2) (for conduction along the y— direction). At this point we 
assume local thermal equilibrium (LTE) which we prove from our simulation data at the end of 
this section. Assuming LTE we write the following approximate form for the energy change AK^ 
during a collision: 

where we have denoted x\°'~'^^ — and AT — Ti — T2. The temperature gradient has been 
assumed to be small and constant. Further we assume that in the close packed limit that we are 
considering, only nearest neighbor pairs {< ij >} contribute to the current in Eq. (6.13) and 
that they contribute equally. We then get the following approximate form for the total current: 

« '-^f . (6.14) 

where Tc is the average time between successive collisions between two particles while is the 
mean square separation along the y— axis of the colliding particles. Finally, denoting the density 
of particles hy p — N/Awe get for the current density: 

^ <i> « . (6.15) 

For strains e^x and eyy in the x and y directions we have p = Po/[(l + ^xx) (1 + ^yy)]- We 
estimate y^ and from a simple equilibrium free-volume theory, known as fixed neighbour free 
volume theory (FNFVT) [ See Sec. 5. 3 of chapter-5 ]. In this picture we think of a single disk 
moving in a fixed cage formed by taking the average positions of its six nearest neighbor disks [see 
Fig. 5. 10]. For different values of the strains we then evaluate the average values [yc]fv and [rdfv 
for the moving particle from FNFVT. We assume that the position of the center of the moving 
disk Po{x,y), at the time of collision with any one of the six fixed disks, is uniformly distributed 
on the boundary B of the free-volume. Hence [yc]fv is easily calculated using the expression: 

[yclfv = -j^ , (6.16) 
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Figure 6.8: Plots showing comparision of the analytically calculated values of [y1\fv and [r^f^, with 
those obtained from a free volume simulation of a single disk moving within the free volume cage. The 
free volume corresponds to a starting unstrained triangular lattice at r/ = 0.85 which is then strained 
along X or y directions. 

where Bi is the part of the boundary B of the free volume when the middle disk is in contact with 
the i^^ fixed disk, ds is the infintesimal length element on B while Lq is the total length of B. Let 
the unstrained lattice parameters be a°, aPy = ■\/3a°/2. Under strain we have = a°(l + e^x) 
and tty = a°(l + tyy). Using elementary geometry we can then evaluate [yl]fv ffom Eq. (6.16) 
in terms of e^xi ^yy and the unstrained lattice parameter a°. An exact calculation of [rjj^ 
is nontrivial. However we expect [rj/t, = c Vjl"^ /T^^"^ where Vf^ is the "free volume" [see 
Fig. 5.10] and c is a constant factor of 0(1) which we will use as a fitting parameter. The 
calculated values for [y1]fv and [rj/t, are shown in Fig. 6.8. Also shown are their values obtained 
from an equilibrium simulation of a single disk moving inside the free volume cage. Thus we 
obtain the following estimate for the heat current: 

c Vfi 

We plot in Fig. (6.3,6.5) the above estimate of [jy]fv along with the results from simulations. 
We find that the overall features of the simulation are reproduced with c = 0.42. 

For small strain [yl{exx)\fv ~ 0.5 + ae^x - Piel^, [yl{(^yy)\fv ~ 0.5 - aeyy + /^se^^ and 
[Tc]fv ~ (71 + 72e - 73e^) where e stands for either exx or eyy and a, (3i, P2, 7i,72,73 are all 
positive constants that depend only on (we do not write them explicitly since the expressions 
are messy and unilluminating). For r] = 0.85 we have a = 7.62, Pi = 121.77, P2 = 124.37 
and 7i = 0.02, 72 = 0.33, 73 = 1.125. From these small strain scaling forms we find that 
jy{eyy) always decreases with positive eyy and increases with negative or compressive eyy (note 
that we always consider starting configurations of a triangular solid of any density). On the other 
hand the sign of the change in jy{exx) will depend on the relative magnitudes of a. Pi and 7 . 
For starting density 77 = 0.85, jy{exx) decreases both for positive and negative exx- 'n Fig. 6.9 
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Figure 6.9: Plot showing effect on jy of negative strains applied in the x and y directions. The 
system is prepared initially in a triangular lattice at r/ = 0.85. Note that negative exx reduces jy whereas 
negative eyy increases i,. 




Figure 6.10: Comparison of simulation results for jy with the approximate formula in Eq. (6.15) 
where Tc and are also calculated directly from the same simulation. The results are for a 40 x 10 
system with starting value of = 0.85 and strained along x— direction. 



we show the effect of compressive strains txx and tyy on the heat current jy and compare the 
simulation results with the free volume theory. 

It is possible to calculate y"^ and Tc directly from our nonequilibrium collision time dynamics 
simulation. The mean collision time Tc is obtained by dividing the total simulation time by the 
total number of collisions per colliding pair. Similarly y"^ is evaluated at every collision and we 
then obtain its average. Inserting these values of and {y'^) into the right hand side of Eq. (6.15) 
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Figure 6.11: Plot of temperature profile and fourth moment of velocity for a strained 40 x 10 
lattice. The unstrained packing fraction was t] = 0.85 and the system was strained to e^x = 0.0625. 



we get an estimate of the current as given by our theory (without making use of free-volume 
theory). In Fig. 6.10 we compare this value of the current jy, for strain e = e^x, and compare it 
with the simulation results. The excellent agreement between the two indicates that our simple 
theory is quite accurate. 

We have also tested the assumptions of a linear temperature profile and the assumption of 
local thermal equilibrium (LTE) that we have used in our theory. In our simulations the local 
temperature is defined from the local kinetic energy density, i.e. ksT = {mu^/2). Local thermal 
equilibrium requires a close to Gaussian distribution of the local velocity with a width given 
by the same temperarture. The assumption of LTE can thus be tested by looking at higher 
moments of the velocity, evaluated locally. Thus we should have (u^) = 8{kBT/mY . From our 
simulation we find out {u^{y)) and ksTiy) as functions of the distance y from the cold to hot 
reservoir. The plot in Fig. 6.11 shows that the temperature profile is appoximately linear and 
LTE is approximately valid. We use our theory only in the solid phase and in this case there is 
not much variation in the direction transverse to heat flow (s— direction). 

Finally we find that the heat conduction in the small confined lattice under small strains shows 
a linear response behaviour. This can be seen in Fig. 6.12 where we plot j,^ versus AT = T2— Ti 
for a 40 X 10 triangular lattice at 77 = 0.85. Note that, as mentioned in the introduction, the 
bulk thermal conductivity of a two-dimensional system is expected to be divergent and the linear 
response behaviour observed here is only relevant for a finite system in certain regimes (solid 
under small strains). 
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Figure 6.12: Plot of jy versus AT = T2 - Ti for a 40 x 10 triangular lattice at r/ = 0.85. We see 
that the current increases linearly with the applied gradient. 

6.3 Summary and Conclusion 

In this chapter we have studied heat conduction in a two-dimensional solid formed from hard disks 
confined in a narrow structureless channel. The channel has a small width (~ 10 particle layers) 
and is long (~ 100 particles). Thus our system is in the nanoscale regime. We have shown that 
structural changes that occur when this solid is strained can lead to sudden jumps in the heat 
current. From the system sizes that we have studied it is not possible to conclude that these 
jumps will persist in the limit that the channel length becomes infinite. However the finite size 
results are interesting and relevant since real nano-sized solids are small. We have also proposed 
a free volume theory type calculation of the heat current. While being heuristic it gives correct 
order of magnitude estimates and also reproduces qualitative trends in the current-strain graph. 
This simple approach should be useful in calculating the heat conductivity of a hard sphere solid 
in the high density limit. 

The property of large change of heat current could be utilized to make a system perform as a 
mechanically controlled switch of heat current. Similar results are also expected for the electrical 
conductance and this is shown to be true at least following one protocol of straining in Ref.[167]. 
From this point of view it seems worthwhile to perform similar studies on transport in confined 
nano-systems in three dimensions and also with different interparticle interactions. 
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